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ABSTRACT 


This  thesis  presents  a  new  approach  to  the  digital  solution 
of  the  power  system  load-flow  problem. 

The  most  popular  approach  in  use  today  is  the  Newton-Raphson 
method  which  is  fast  but  uses  a  considerable  amount  of  core 
memory . 

The  method  proposed  here  is  a  modification  of  the  Newton- 
Raphson  method,  which,  although  slower  in  convergence,  can  save 
considerable  time  in  solution  as  it  eliminates  the  need  for  re¬ 
calculation  of  the  coefficients  of  the  linearized  system  at  each 
iteration. 

The  core  requirements  are  only  slightly  higher  than  Newton- 


Raphson'  s  . 
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I.  THE  LOAD-FLOW  PROBLEM 


1.1  DEFINITION 

The  load  flow  problem  is  the  problem  of  determining  the 
voltages  and  power  flows  in  a  power  system  for  a  given  set  of  loads 
and  system  connections.  It  is  necessary  to  know  in  advance,  the 
voltages  and  generator  outputs  required  to  carry  the  loads  in  a 
satisfactory  manner. 

The  terminal  or  bus  conditions  in  a  network  of  distribution 
are  likely  to  change  continuously.  Indeed,  the  demand  of  power  by 
the  users  is  variable,  and  the  generation  has  to  be  changed  accord¬ 
ingly.  Different  patterns  of  load  at  the  load  buses,  require 
different  patterns  of  generation  at  the  generation  buses.  However, 
not  all  patterns  of  generation  are  possible  because  of  limitations 
in  bus  voltages  and  line  current  carrying  capacity. 

Each  pattern  of  generation,  for  a  given  pattern  of 
loads,  implies  a  cost  of  generation  plus  a  certain  amount  of  power 
losses  in  the  network.  If  an  optimum  economical  operation  is  to 
be  implemented,  several  patterns  of  generation  have  to  be  tried,  so 
that  the  best  pattern  can  be  found- 

Similarly,  in  order  to  ensure  that  the  powers  carried  and  the 
voltages  are  within  the  safety  levels,  a  number  of  studies  have  to 
be  made  with  different  generation  patterns.  This  also  must  be  done 
for  testing  the  electric  stability  of  the  power  system  when  operated 


with  a  specified  pattern  of  loads  and  generations. 

Since  the  power  demand  is  increasing,  new  power  stations  and 
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new  lines  of  interconnection  have  to  be  built  up.  Preliminary 
planning  studies  have  to  be  carried  out  to  check  the  feasibility 
of  the  new  supply  as  far  as  stability,  safety  levels  and  economic 
optimality  are  concerned.  The  performance  of  a  power  system  is 
known  completely  if  we  know  its  load-flow  solution,  and  from  this 
stems  the  importance  of  the  load-flow  studies. 

The  essential  problem  is  that  of  determining  the  node  voltages 
in  the  system  as,  once  these  are  known,  the  power  flows  may  be 
calculated  by  use  of  equation  (1.1), 


P  -jQ  =  E  *(E  -  E  )y  +  E  *E  — 23. 

pq  pq  p  p  q  pq  p  p  2 


d.i) 


where  p  and  q  stand  for  the  code  number  given  to  two  of  the  buses 
of  the  network.  The  current  flowing  into  the  system  at  bus  "p", 
can  be  calculated  by  equation  (1.2) 


I 

P 


(1.2) 


A  single  phase  representation  is  enough  because  all  problems 
dealt  with  are  symmetric  and  balanced. 

Several  representations  of  the  network  are  possible.  Each 
representation  has  several  associated  methods  of  solution  of  the 
load  flow  problem. 

In  1929  the  alternating  current  network  analyzer,  which  is 
a  special  purpose  analogue  computer,  was  created.  This  device  was 
used,  among  other  studies,  for  the  solution  of  load-flow  problems. 


. 
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The  network  analyzer  severely  limits  the  size  of  problem 
which  may  be  handled,  requiring  a  simplified  representation  of  the 
system,  in  many  cases.  However  until  the  mid  1950s  digital  computers 
were  not  satisfactory  due  to  core  limitations.  From  the  late  1950s 
onward,  the  progress  in  digital  computer  performances  has  kept 
pace  with  the  increasing  complexity  and  size  of  power  system 
problems.  Nowadays  digital  computers  have  replaced  alternating 
current  network  analyzers  as  the  tool  for  the  power  system  engineer. 
At  the  same  time  the  amount  of  memory  required  and  the  computer  time 
spent  have  become  the  indications  of  goodness  of  a  load  flow  method. 
This  is  because  the  load-flow  problem,  as  it  has  been  shown  before, 
has  to  be  solved  many  times,  and  both  the  amount  of  computer  memory 
and  time  used  affect  costs. 

1.2  METHODS  OF  SOLUTION 

Different  mathematical  descriptions  of  the  network 
result  in  different  load-flow  equations.  The  first  method  used 
for  digital  solution  implemented  the  loop  frame  of  reference  in 
admittance  form.  However,  this  form  takes  a  longer  data  preparation 
than  the  ones  that  use  the  bus  frame  of  reference,  in  either 
admittance  or  impedance  form. 

Currents  or  voltages  can  be  chosen  as  independent  variables. 

If  currents  are  selected,  the  impedance  form  must  be  used.  Equation 
(1.2)  along  with  (1.3)  are  used 


(1.3) 


_  1+  _ 


R 

where  E  p  is  the  voltage  at  the  reference  bus. 

If  the  admittance  form  is  selected,  the  equations  used  have 
the  voltages  as  independent  variables,  and  several  equations  such 
as  (1.4)  and  (1.5)  could  be  used  depending  on  the  method  of  solution 
applied 


Ep 


=  J_  (gp.  -  jQp 

Ypp  v  Ep* 


-  ^pqEq) 

q^p 


p  =  1,...,N  ,  p  ^  s  ("s"  means 

slack  bus,  see  below)  (1.4) 


Pp  =  |{ep(epGpq  +  fqBpq)  +  fp(fqGpq  -  eqBpq) } 


Qp  =  2{fp(eqGpq  +  f qBpq)  -  ep(fqGpq  -  eqBpq)}  (1.5) 

where  (Ypq  =  Gpq  -  jBpq) 

There  are  four  quantities  associated  with  each  bus.  These  are 
the  real  and  reactive  power,  and  the  real  and  reactive  voltage.  Only 
two  of  the  quantities  are  known  at  each  bus  before  the  solution  of 
the  problem.  Depending  on  which  quantities  are  known,  the  bus  is 
called  differently  and  processed  in  the  appropiate  way.  We  can  have 
the  following  kinds  of  buses. 

Power  buses 

-  generation  buses 

-  load  buses 

-  generation  plus  load  buses 

Voltage  buses 

-  voltage  controlled  buses 

-  slack  buses 

Power  buses  are  the  ones  at  which  the  powers  flowing  into 
the  system  or  out  of  the  system  are  known,  and  the  active  and 


-  5  - 


and  reactive  voltage  are  unknown.  If  the  power  at  the  bus  is  flow¬ 
ing  into  the  system,  we  call  it  a  generation  bus.  If  the  power  is 
flowing  out  of  the  system  we  call  it  a  load  bus.  When  there  is 
generation  and  load  at  the  same  bus,  it  is  called  generation  plus 
load  bus;  in  general  the  generation  and  the  load  in  such  buses  are 
not  balanced.  Therefore  at  such  buses  a  certain  amount  of  power 
is  flowing  into  the  system  or  out  of  it.  The  powers  have  a  positive 
sign  when  flowing  in,  and  negative  when  flowing  out. 

Voltage  buses  are  the  ones  where  either  the  voltage  magnitude 
is  constant  or  both  components,  the  real  and  the  reactive  voltage 
are  constant.  The  buses  at  which  the  voltage  magnitude  is  constant 
can  change  their  active  and  reactive  voltages  provided  the  voltage 
magnitude  is  kept  constant.  In  such  buses  the  real  power  is  known 
and  constant,  and  the  reactive  voltage  varies  so  that  the  voltage 
magnitude  is  constant.  These  are  called  voltage  controlled  buses. 
When  both  components  of  the  voltage  are  known,  the  bus  is  called 
a  slack  bus.  In  that  case  neither  the  real  power  nor  the  re¬ 
active  power  are  known  before  the  solution  of  the  load-flow  problem. 

Trial  and  error  techniques  can  be  used  to  account  for  such 
special  cases  as  limited  available  reactive  power  at  voltage 
controlled  buses,  tap-changing  under  load  transformers  and 
specified  net  interchange  of  power.  Although  it  takes  a  longer 
time,  the  solution  is  always  reached. 

The  losses  in  a  power  system  are  a  function  of  the  node  or 
bus  voltages,  which  are  unknown  before  the  load-flow  solution. 


V 
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Since  the  generations  have  to  match  the  loads  plus  the  losses, 
there  must  be  at  least  one  slack  bus,  with  an  indeterminate  power 

inflow,  in  order  to  generate,  at  that  node,  as  much  power  as  need¬ 
ed  to  balance  the  losses. 

A  very  simple  example  is  given  now  to  show  the  different 
types  of  buses  and  the  equations  to  be  solved.  The  admittance  form 
and  equation  (5)  are  chosen  for  the  mathematical  formulation.  The 
problem  consists  of  four  buses  as  shown  in  figure  1-1. 


Figure  1.1  SIMPLE  EXAMPLE 

Bus  1  is  a  slack  bus.  Its  voltages  eq,  fq  are  fixed  and  known. 
Its  powers  Pq,  Qq,  are  unknown  and  will  be  computed  at  the  end  of  the 
process,  once  the  voltages  at  all  other  nodes  are  known,  by  means  of 
equation  (1.5). 

Bus  2  is  a  generation  plus  load  bus.  Its  voltages  ,  ^2  are 
unknown  and  are  found  by  solution  of  the  load-flow  problem.  The 
generation  at  this  node  is  P2  “  j Q2  anc*  its  load  L2  -  1^2  *  The 


’ 
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power  flow  into  the  system  is  thus:  (P2  -  L2)  -  j  (Q2  “  M2).  The 
way  in  which  equation  (1.5)  is  applied  in  this  case  is  the  following, 

P2  -  L2  =  e2(e1G21  +  f l®2l)  +  f2^1^21  "  el®2l)  + 

+  e2^e3^23  ^3®23^  ^2^3^23  —  e3^23^  (1*6) 

(the  symbols  with  a  dash  are  known  quantities) 

Q2  -  M2  =  f 2  (elG21  +  f  1^2l)  "  e2^1^21  "  el®21^  + 

+  f 2  (£3623  +  f 3^23^  ”  e2  ^ 3*^23  ~  e3^23^  (1*7) 

Bus  3  is  a  load  node.  The  power  outflow  or  negative  in  - 
flow  is  -(L3  -  jM2) 

-  L3  =  e3(e1G31  +  fiB31)  +  f 3 (f  1G31  -  eiB31)  + 

*4*  63 (e2G32  "4*  £3^32^  "£  ^3^^2^32  —  ^2^32^  "£ 

+  63 (e^G3^  +  £4^34)  +  f3(f4G34  -  e^i^)  (1.8) 

-  M3  =  f3(e1G31  +  fiB31)  -  e3(f1G31  -  exB31)  + 

+  f3(e2G32  +  £2^32)  “  e3 2^32  ”  e2®32)  + 

+  f 3 (e4G34  +  £4^34)  ~  e3 (f 4^34  “  e4®34) 


Bus  4  is  a  voltage  controlled  bus.  The  voltage  magnitude  is 


kept  at  the  value  V4,  hence, 


■ 


■*  I 
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V42  =  ^42  +  f42 


(1.9) 


the  reactive  power  needed  will  be  calculated  once  the  solution  has 
been  reached  by  means  of  equation  (1.5).  As  for  the  active  power, 
we  do  not  have  in  this  case  any  power  inflow  at  that  node,  since 
there  is  no  real  power  generation  nor  real  load,  therefore, 


0  -  e^Ce^G^  +  fiB^)  +  f^f]/^  -  6164!)  + 


+  e4 (e3G43  +  f3B43)  +  f4^f3G43  “  e3B43^ 


(1.10) 


These  six  equations  have  to  be  solved  to  find  e2 ,  f2>  63, 
f3>  ^4>  ^4*  As  these  equations  are  nonlinear  equations,  an 
iterative  technique  must  be  used.  Several  methods  could  be  used 
for  solving  this  set  of  nonlinear  equations  by  means  of  a  digital 
computer . 

1.2.1  CLASSIFICATION  OF  THE  METHODS  OF  SOLUTION 

As  a  classification  of  the  methods  of  solution  of  the  load- 
flow  problem,  that  given  by  Laughton  and  Humphrey  Davies  (1)  of  direct 
load  flow  methods  and  iterative  load-flow  methods,  could  be  taken. 

Sasson  and  Jaimes  (2),  came  up  later  with  a  new  classification, 
(nodal  admittance  matrix  methods,  variational  matrix  methods,  and 
nodal  impedance  matrix  methods),  objecting  to  the  ambiguity  of  the 
former  classification  as  the  direct  methods  are  iterative  in 
nature. 

This  objection  could  be  refuted  by  changing  the  name  "iterative 


methods".  Since  all  of  the  so  called  iterative  methods  correct 
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the  voltage  of  only  one  bus  at  a  time,  they  could  be  called  "bus 
by  bus"  methods,  and  this  name  will  be  used  in  the  rest  of  this 
section. 

In  a  recent  compendium  of  computer  methods  in  power  system 
analysis,  Stagg  and  El-Abiad  (3) ,  describe  the  most  usual  load- 
flow  methods  without  any  classification.  However  ,  the  methods 
presented  can  all  be  fitted  into  the  classification,  of  direct 
methods  -  bus  by  bus  methods. 

bus  by  bus  methods  Gauss  iterative  method  using  Y^us 

Gauss  iterative  method  using  Zbus 
Gauss-Seidel  method  using  Y^us 
Gauss-Seidel  method  using  Zbus 
Relaxation  method  using  Y^us 

direct  methods  Newton-Raphson  method  using  Ybus 

Gauss  iterative  method  using  Y^OQp 

A  unified  nonlinear  programming  approach  to  solve  the  load- 
flow,  minimum  loss,  and  economic  dispatch  problems  has  been  shown 
to  be  possible  by  Albert  M.  Sasson  (4) . 

All  methods  start  with  a  first  assumption  of  the  bus  voltages. 
Bus  by  bus  methods  use  an  algorithm  which  then  computes  a  better 
approximation  for  each  bus  on  a  single  bus  basis.  The  order  of 
computation  follows  the  given  numbers  of  the  buses  in  Gauss-Seidel 
methods  and  Gauss  iterative  methods.  In  the  relaxation  method, 
the  only  bus  whose  voltage  is  corrected  at  each  step,  is  that  with 
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the  largest  residual,  which  is  a  measure  of  the  mismatch  in  power. 

In  the  direct  methods  the  correction  of  the  voltages  come 
from  an  operation  which  uses  all  the  available  information.  This 
operation  is  the  solution  of  a  linear  system  of  equations  in  the 
Newton-Raphson  method,  and  is  the  calculation  of  several  vectors 
and  matrix  multiplications  in  the  case  of  the  Gauss  iterative 
method  using  Yloop. 

1.2.2  FEATURES  OF  THE  MOST  IMPORTANT  METHODS 

There  is  not  a  best  method  in  load-flow  because  there  can 
always  exist  a  specified  power  system  for  which  one  of  the  methods, 
which  does  not  work  as  well  as  the  other  methods  with  other 
problems,  gives  better  results  than  the  others.  Moreover,  even  by 
changing  the  initial  guess  in  the  same  lead-flow  problem,  it  is 
possible  to  find  different  performances  of  the  methods.  All  methods 
have  advantages  and  disadvantages  and,  in  most  of  the  cases,  the 
selection  of  one  method  or  another  depends  on  the  weight  placed 
on  every  advantage  and  disadvantage. 

Bus  by  bus  methods  in  general  have  the  disadvantage  that  in 
order  to  have  a  small  mismatch  in  power,  they  have  to  iterate  until 
a  considerably  smaller  increment  in  voltage  is  reached. 

The  early  method  of  Ward  and  Hale  (5) ,  has  been  improved  by 
other  more  sophisticated  bus  by  bus  methods.  In  their  book  (3), 
Stagg  and  El-Abiad,  (p.  328),  show  how  the  Gauss-Seidel  method  using 
Ybus  takes  the  lowest  time  per  iteration.  But  as  can  be  seen 
later,  (pp.  330  and  331),  this  advantage  is  balanced  by  the  high 
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number  of  iterations  for  solution. 

As  far  as  time  is  concerned,  it  can  be  seen  that  Newton-Raphson 
method  yields  the  best  results,  but  there  are  other  factors  to  account 
for.  One  of  these  is  the  computing  time  for  both  processing  system 
data  in  order  to  obtain  the  parameters  for  the  iterative  solution, 
and  modifying  network  data  and  introducing  operating  changes.  In 
this  regard  the  methods  which  use  the  bus  admittance  form  are  better 
off  than  any  other  form,  because,  provided  that  mutual  coupling  is 
not  considered,  the  computation  and  changes  of  Y|3US,  are  very  easy 
and  straightforward. 

There  is  still  another  factor,  the  storage  required.  The 
Gauss-Seidel  method  using  Y^us  needs  less  computer  memory  than 
Newton-Raphson  method  using  Ybus.  Whether  the  extra  storage  is 
offset  by  the  savings  of  time  by  the  Newton-Raphson  method,  has  to 
be  decided  by  the  utility  companies  which  use  such  methods. 

Tinney  and  Hart,  (6),  have  obtained  some  improvements  in  the 
practical  application  of  the  Newton-Raphson  method  to  real  power 
systems.  Their  work  along  with  other  studies,  such  as  the  work 
by  Feingold  and  Spohn  by  Electricite  de  France,  (7),  seems  to 
indicate  that  Newton-Raphson  is  the  method  most  commonly  employed 
by  the  big  utility  companies. 

1.3  BRIEF  SCOPE  OF  THE  THESIS 

An  improvement  in  time  with  respect  to  the  Newton-Raphson 
method  has  been  found. 

A  different  approach  to  the  solution  of  the  set  of  nonlinear 
equations  (1.5),  leads  to  a  different  iterative  method  of  solution. 
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The  main  features  of  this  method  are: 

a)  The  calculation,  once  only,  of  the  non-diagonal  coefficients 
of  the  linearized  set  of  equations.  With  the  Newton-Raphson 
method  this  calculation  must  be  made  for  each  iteration. 

b)  The  feasibility  of  the  implementation  of  fast  approximate  so¬ 
lution  of  the  linearized  system. 

One  of  the  advantages  over  the  Newton-Raphson  method  is  that 
not  much  time  is  spent  during  the  early  stages  for  the  solution 
of  the  linearized  system,  since  it  is  going  to  yield  results, 
which  are  still  far  from  the  true  solution.  An  exact  time  con¬ 
suming  method  is  not  implemented  until  the  voltages  are  close  to 
the  final  solution.  This  idea  had  been  already  suggested  by  Van 
Ness  in  1959,  (8),  but  it  had  not  been  developed. 


. 


II. 


THE  NEWTON-RAPHSON  METHOD 


2.1  DESCRIPTION 

As  the  Newton-Raphson  method  is  presently  the  most  popular 
method  in  use,  it  will  be  described  in  some  detail.  Alternative 
methods  may  then  be  compared,  with  regard  to  speed  and  accuracy, 
with  the  Newton-Raphson  method. 

The  Newton-Raphson  iterative  technique  for  finding  the  root 
of  a  function  is  used  by  the  Newton-Raphson  method  for  solving  the 
set  of  equations  (1.5).  This  method,  as  it  is  known  to  the  power 
systems  engineers,  is  nothing  but  the  generalization  of  the  Newton- 
Raphson  technique  for  finding  roots  of  a  one  dimensional  function, 
to  a  multidimensional  space. 

The  set  of  equations  (1.5),  may  be  written  in  a  more  general 
fashion  (2.1) 

pp  ~  §p(®l»®2>,,,»  eN >  f 1 >  f 2  > • • • ’ fN) 

Qp  —  hp  (eq,  e2  >  •  •  •  >  sjq  > p  1  >  2  ’  *  ’  *  *  ^N^  (2.1) 

p  =  1 , . .  .  ,N 

If  the  set  of  functions  gp  and  hp  are  expanded  about  the  point 
(e°q,  e°2 , • • • ,e°N, f°q ,f°2 , . . . , f °n)  in  a  2xN  dimensional  coordinate 
system,  in  a  Taylor  series,  and  drop  the  second  order  derivatives 
and  those  of  higher  order,  then 

pp  =  gp(e°l,e02> . . • ,e°N,f°i,f°2, • • • ,f°N>  + 


' 


. 


-  Ik  - 


+  -Sp 

6ep 


(ei  -e°^)  +  ...  +  |||| 
o  ^  'o 


(ei-e0^)  +  ...  +  ||®j 


+  -Sp  s  (f-i  -f°  )  +  +  -§P  i  +  +  ^§E  j 

6fi|0Ul  £  1;  +  <$f  i  jo  ^ 1  1  +  6fN|o 


Qp  =  hp(e°1,e02, • . . ,e0N,f01}f02, • • . ,f°N)  + 


+  — E  j  fp.-pO,)  +  +  ^!}e|  (e  -e° . )  +  +  --P 

6elloC  1  1;  *”  <5ei|0  {ei  e  +  •••  5eN 


+  ~e  (fi-fv  +  ...  h-  §1*1  +  ...  + 


6f  ■ 


Nb 


Calling 


i  ^  s,  (S  slack  buses,  see  App. 


P°  =  gp(e°1,e02, • • • ,e°N,f01,f°2, . . . ,f°N) 


Q°p  =  hp(e°1,e°2, . . . ,e0N,f01}f°2, • • • ,f°N) 


iO. I  „  ^SE 
3  P1  6e± 


,o  .II 


<5g. 


Pi  6f. 


ao  .III  =  ^hpi 
J  P1  6ei|0 


ao  .IV  „  ^hE| 

J  P1  &f±l 


where  the  subscript  "o"  means  at  the  point  (e0-^, .  . .  ,e°N, 


(eN-e°N)  + 

(fN-f°N)  +  £ 

(eN-e°N)  + 

(fN-fV  +  £ 

6) 


f°l,...f°N) 


Equation  (2.2)  can  be  rewritten  in  matrix  form 
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•  •  • 

»  •  •  •  •  •  •  •  •  j  »  «  .  •  •  e  o*o 

.  a  . 

Pp  -  p°p 

•  •  a 

i°  .!  ,  -?o  I 

J  pi  ...  |  ...  J  pi 

ei  “  e°i 

= 

•  oil  •  i  •  |  iii  §  o  •  ©  .  © 

— 

•  .  « 

Qp  -  Q°P 

•  •  •  •  •  •  a  •  •  |  ©  •  •  •  •  •  ©  •  © 

iO  .III  .  AO  .IV 

...  j  ...  ...  j  pi  .  .  . 

... 

f  .  _  f  o  . 

r  l 

•  •  • 

—  _ 

•  •  •  ii*  ...  |  ...  ...  ... 

»  .  * 

p  =  1, . . . ,N  ;  p  ^  s  (2.3) 

i  =  1,...,N  ;  i  ^  s  ,  (S  slack  buses) 


It  must  be  pointed  out  that  equation  (2.3)  is  an  approximate 
equation  because  the  second  order  derivatives  and  those  of  higher 
order,  have  been  neglected.  It  can  be  noticed  that  no  equations 
for  the  slack  buses  are  taken,  (see  appendix  6  for  more  than  one 
slack  bus) . 

As  a  result  of  the  first  solution  of  equation  (2.3)  an  incre¬ 
ment  of  voltage  is  obtained,  which  will  bring  the  voltages  to  their 
new  value.  Indicating  the  first  increment  by  A^ei  and  A^fi 

Al  e-jL  —  0  i  —  e°i  —  ^^~i  —  ^o  ^ 

Aifi  =  f±  -  f°i  =  f1!  -  f°i  (2.4) 

The  superscript  "1"  is  used  to  indicate  the  voltages  after 
the  first  solution  of  (2.3),  or  first  iteration.  Of  course 

e-^i  =  e°i  -1-  A^-e^ 


£-q  =  f°i  +  Aifi 


(2.5) 
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It  can  be  said  that,  in  a  general  iteration,  the  superscript 
of  the  increment  in  voltage  is  the  number  of  the  iteration  whereas 
the  superscript  of  the  increment  in  power  is  that  of  the  former 
iteration.  The  power  increments  are  calculated  with  the  voltages 
obtained  in  the  iteration  before,  and  the  same  for  the  terms  of 
the  jacobian,  jkpj^,  etc.,  where  "k"  indicates  that  these 

terms  are  calculated  with  the  voltages  eki,  f^i,  obtained  at  the 
"kth"  iteration.  The  general  form  of  equation  (2c3),  that  which 
would  yield  the  "k+1"  set  of  voltages,  at  the  "k+l"  iteration,  is 
the  next 


AkP 

jk  I  ]  jk  II- 

Ak+1e 

AkQ 

jk  III  l  jk  IV 

1 

X 

fik+lf 

The  same  process  is  followed  at  each  iteration.  When  a 
new  set  of  voltages  e^p5fkp  are  obtained,  the  functions  gp  and 
hp  (2.1)  are  expanded,  in  a  Taylor  series  about  the  point  (ek^, 
ekg  , .  .  .  ,  ek^ ,  f^i , f^g  ,  •  .  •  , fS\j)  >  and  then  proceed  to  the  next  step 
as  in  the  first  iteration. 

AkP  and  AkQ  become  smaller  and  smaller  as  the  number  of 
iterations,  "k",  increase.  When  these  increments  become  less 
than  a  specified  tolerance,  the  process  may  be  stopped  and  e^-j_, 
fki>  considered  as  the  solution.  A  flow  chart  of  the  method  can 
be  found  in  appendix  3. 

It  could  be  easily  shown  that, 


' 


-  IT  - 

jIik  =  eiGik  "  fiB±k  i  ^  k 
G^ik  =  e±Bik  +  fiGik  i  ^  k 
Jinik  -  JIXik 
JIVik  =  -  JXik 

J^ii  =  epGpp  -  +  ^(ekGj_k  +  fkBik) 

J^ii  =  eiBii  +  fiGii  +  ^(^kGik  ”  akBik) 

J^^ii  =  eiBii  +  fiGn  -  ^(fk^ik  “  ®kBik) 

JIVii  =  fiBii  ~  eiG^j_  +  ^(ekGik  +  fkBik) 

so  the  linearized  system  to  be  solved  would  be  (2.7), 

APp  =  5^eq(epGpq  “  ^pBpq)  +  Afq(epBpq  +  fpGpq)  }  + 

+  AepE(eqGpq  +  fqBpq)  +  AfpZ(fqGpq  -  eqBpq) 

AQp  =  <|^eq(epBpq  +  fpGpq)  +  Afq(fpBpq  -  epGpq)  }  + 

+  AepZ(eqBpq  -  fqGpq)  +  AfpL(eqGpq  4-  fqBpq) 

p  f  1,2,. ..,N  (p  *  S)  (2.7) 

Section  5.5  of  P.  Henrici’s  book  on  numerical  methods  (10), 
gives  a  fair  insight  into  the  solution  of  a  set  of  nonlinear 
equations  by  the  Newton-Raphson  method.  For  further  details, 
section  3.3  of  Isaacson  and  Kellers’  book  (11)  can  be  consulted. 

2.2  NUMERICAL  PROPERTIES 


It  is  proven  in  section  3.3.2  of  reference  (11),  that  if 
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the  following  set  of  conditions  are  satisfied,  the  Newton- 
Raphson  method  converges  to  a  solution. 

Calling  J  the  jacobian  and  E°,F°  the  initial  iterate  of 
the  voltage,  it  can  be  proven  that, 


-  if  J~1(E°,F°)  exists  and  its  norm  is  bounded,  "a"  being 
the  bound,  and 

-  if  the  norm  of  the  first  set  of  increments  A^ei,  A-'-fi 
has  "b"  as  a  bound  for  its  norm,  and 

-  if  for  all  the  E,F  that  lie  in  the  2b-sphere  about  E°,  F° 
any  of  the  functions  (2.1)  gp(E,F)  or  h  (E,F)  have  con¬ 
tinuous  second  order  derivatives  with  respect  to  any 


pair  of  components  of  E,F,  which  satisfy: 
%r,(E.F'),  ^  c 

and 


LbvL EjH  [  < 

<5ej  6ek  1  2xN 


_  if 


axbxc  -  — 


then , 


all  Newton-Raphson  iterates  are  within  the  2b-sphere  about  E°,F°, 
and  these  iterates  converge  towards  a  vector  which  satisfies  the 
equations  (1.5). 

Furthermore,  if  the  function  (2.1)  gp  and  hp  have  two 
derivatives  and  the  jacobian  is  non-singular  at  the  solution, 
the  convergence  of  Newton-Raphson * s  method  is  quadratic,  which 
means  that  the  error  at  the  "k+1"  iteration  is  less  than  the 
square  of  the  error  at  the  "k"  iteration  multiplied  by  a  bounded 


constant . 
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2.3  ANALYSIS  OF  COMPUTATION  TIME 

For  analysis  of  the  performance  of  various  methods  a  number 
of  test  systems  have  been  chosen.  The  IEEE  14,  30,  and  57  bus 
test  systems,  along  with  other  test  systems  with  higher  number 
of  buses,  were  put  forward  by  the  American  Electric  Power 
Conference  for  investigators  to  use  in  their  research  concerned 
with  electrical  power  systems. 

These  test  systems  have  the  average  characteristics  of 
existing  power  systems  as  far  as  electric  stability  and  numerical 
properties  are  concerned. 

The  data  for  the  IEEE  14,  30,  and  57  bus  test  systems  can 
be  found  in  reference  (9).  Applying  the  Newton-Raphson  method 
to  these  three  cases,  it  is  found  that  3  complete  iterations  and 
the  start  of  a  fourth  iteration  are  needed. 

The  iterations  have  been  split  into  three  parts  as  far  as 
time  measuring  is  concerned.  These  subdivisions  are: 

-  Calculation  of  new  currents  and  powers,  and  checking  of 
the  differences  between  the  computed  powers  according  to 
the  iterated  voltages,  and  the  scheduled  powers. 

-  Calculation  of  the  new  terms  of  the  Jacobian  and  the  new 
set  of  differences  in  power. 

-  Solution  of  the  linear  system  and  correction  of  the 
voltages  with  the  set  of  differences  in  voltage  obtained 
from  the  linear  system. 

The  times  used  are  those  corresponding  to  the  program  shown 
in  appendix  3.  The  times  are  average  times  because  the  length 
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Figure  2 . 2 


TIME  PER  ITERATION  VS.  NUMBER  OF  BUSES 


-  21 


of  the  same  problem  is  not  a  constant,  due  to  the  fact  that  the 
computer  used  utilizes  different  modes.  The  unified  scale 
times  are  the  following: 

TABLE  2.1  TIMES  FOR  STAGES  OF  COMPUTATION  IN  THE  NEWTON- 


RAPHSON  METHOD 

n°  of  buses 

time  per 

new  currents 

new  terms 

linear  sys. 

iteration 

and  powers 

in  jacobian 

and  new  volts. 

14 

0.3254  sec. 

0.0166  sec. 

0.1548  sec. 

0.1540  sec. 

(5.1  %) 

(47.6%) 

(47.3%) 

30 

2.8653  sec. 

0.0526  sec. 

1.5604  sec. 

1.2523  sec. 

(1.8%) 

(54.5%) 

(43.7%) 

57 

19.2078  sec. 

0.1502  sec. 

11.2905  sec. 

7.7670  sec. 

(0.08%) 

(58.7%) 

(41,2%) 

In  brackets  are  the  percentages  with 

respect  to  the 

total 

time  per  iteration.  Graphs 

of  the  times 

and  percentages 

versus 

the  number 

of  buses  are  given  in  figures 

2 . 1  and  2.2. 

The  equation  that  approximately  matches  the  shape  of  the 

2  9 

curve  in  figure  2.2,  is  t  =  KxN  ,  where  t  stands  for  time, 

K  is  a  constant,  and  N  stands  for  the  number  of  buses. 

Most  of  the  time  of  each  iteration,  is  involved  in  the 
calculation  of  the  terms  of  the  jacobian  and  in  the  solution  of 
the  linearized  system.  Any  savings  in  time  by  other  methods  must 
replace  these  parts  of  the  computation  by  more  efficient 
algorithms.  Elimination  of  the  calculation  of  the  jacobian  would 
appear  to  be  the  most  promising  improvement. 


III. 


THE  PROPOSED  METHOD 


3.1  an  ITERATIVE  METHOD  FOR  SOLVING  A  SET  OF  NONLINEAR  QUADRATIC 
EQUATIONS 

H.S.  Wall  (13),  developed  a  modification  of  Newton's  method 
for  approximating  the  roots  of  a  function  f(x). 

xp(i  =  0,1, 2,...) ,  being  the  first  estimation  and  the  succes¬ 
sive  iterates  of  a  root,  the  modification  consists  in  finding  the 
parabola  through  xi,f(xp),  having  the  same  first  and  second  deriva¬ 
tives  as  f(x).  The  intersection  of  the  parabola  through  x-[,f(xj_), 
with  the  x  axis,  is  the  new  iterate  xf+p.  The  Newton-Raphson  me¬ 
thod  finds  new  iterates  by  means  of  the  tangents  to  f(x)  through 
the  points  xp,f(xi). 

The  iterative  technique  implemented  in  this  thesis  finds  a  pa¬ 
rabola  through  x0,f(x0)  having  the  same  first  and  second  derivati¬ 
ve  as  f(x)  at  xQ ,  and  takes  as  a  new  iterate  the  intersection  of 
this  parabola  with  the  x  axis.  The  new  iterates  are  used  for  impro¬ 
ving  the  coefficient  of  the  second  derivative.  As  this  coefficient 
is  unknown,  it  is  taken  as  zero  during  the  first  iteration,  which 
makes  the  first  iteration  equivalent  to  an  iteration  of  the  Newton- 
Raphson  method. 

The  new  method  is  different  from  Wall's  method  in  that  all  the 
parabolas  pass  through  the  initial  estimate  x0,f(x0),  whereas  in 
the  method  by  Wall  the  parabola  at  each  iteration  pass  through  the 
point  xi,f(xi).  Moreover,  the  approximation  used  for  the  coeffi- 
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cient  is  different  for  each  method. 

The  new  iterative  method  is  suitable  for  finding  the  roots  of 
a  quadratic  function, -because  such  functions  have  no  third  deriva¬ 
tive,  and  the  new  method  uses  only  the  first  and  the  second  deri¬ 
vative  . 

The  equations  that  describe  the  load-flow  problem  are  quadra¬ 
tic  with  respect  to  the  bus  voltages,  which  are  the  unknowns. 

Therefore,  a  solution  can  be  expected  from  the  implementation  of 
the  proposed  method  to  the  load-flow  problem. 

The  new  method  has  some  similarities  with  the  Newton-Raphson 
method.  However,  due  to  the  nature  of  the  equations  and  iterations 
of  the  new  method,  its  numerical  properties  are  more  difficult  to 
prove  than  those  of  Newton-Raphson ' s .  The  complete  numerical 
features  of  the  new  method  for  one-dimensional  cases,  are  given  in 
section  3.1.1  •  The  comparison  with  the  one-dimensional  case  for 
the  Newton-Raphson  method,  is  given  in  section  3-1.3  . 

A  sufficient  condition  for  convergence  is  the  multidimensional 
case  is  given  in  section  3-3  •  This  convergence  condition  is  the 
application  of  the  general  convergence  condition  for  the  iterative 
methods  of  solution  of  systems  of  nonlinear  equations. 

In  section  2.2  are  the  requirements  under  which  these  conditions 
are  fulfilled  by  the  Newton-Raphson  method.  Proofs  can  be  found  in 
the  numerical  analysis  references  (10)  and  (ll).  As  for  the  new  method 
in  the  multidimensional  case,  the  requirements  under  which  the  general 
convergence  condition  is  fulfilled,  have  not  been  found  theoretically, 


- 
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because  of  the  complexity  of  the  equations.  However,  the  experience 
with  the  one-dimensional  case,  and  the  arguments  presented  in 
section  3.3  make  one  believe  that,  as  for  the  one-dimensional  case, 
the  requirements  are  the  same  as  for  the  Newton-Raphson  method.  The 
IEEE  lb,  30  and  57  bus  test  systems  show  that,  when  the  proposed 
method  is  used,  the  convergence  follows  the  same  rule  as  the  one¬ 
dimensional  case. 

3.1.1  THE  NEW  ITERATIVE  METHOD  IN  THE  ONE-DIMENSIONAL  CASE.  Con¬ 
sider  a  quadratic  function  q(x),  from  which  roots  are  to  be  found 
by  an  iterative  method.  If  the  value  xQ  is  taken  as  a  first  assump¬ 
tion  of  the  root  x#,  and  q(x)  is  expanded  in  a  complete  Taylor  se¬ 
ries  expansion  about  xQ ,  q(x#)  becomes  (3.1a) 

q(x#)  =  0  =  q(xQ)  +  (x#  -  x0)q'(xQ)  +  ?j(x#  -  x0)2q"(x0) 

(3.1a) 

The  expression  (3- la)  is  exact  since  there  are  no  more 
derivatives  beyond  the  second,  because  q(x)  is  quadratic. 

If  the  nature  of  q(x)  is  known,  there  is  an  easier  way  to  find 
the  exact  development  about  xQ.  Let  us  assume  that  q(x)  is  given  by 

(3.1b), 


q(x#)  =  0  =  (x#  -  rx)  (x#  -  r2)  (3.1b) 

where  r-^  and  r2  are  the  roots  of  q(x)  .  Equation  (3.1b)  can  be 
written  in  an  equivalent  way. 
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qU* 

)  =  qU* 

-  xo 

+ 

X 

o 

u 

(x*  -  XQ  +  xo  - 

i 

* 

X 

i — i 

u 

xo  +  xo  ~  r2> 

(3.1c) 

which  is 

equivalent  to 

(3. Id) 

q(x* 

)  =  (xQ  • 

-  r1)  ( 

xo  "  r2 

+ 

1 — 1 
5m 

1 

o 

X 

+ 

(xo  -  r2 

)HX*  "  xo)  + 

+  (x*  - 

~  xo)2 

(3. Id) 

and 

since 

q(x)  is 

as  in 

(3.1b) 

5 

q(xQ)  = 

(xQ  - 

o 

X 

i — 1 

-  r2) 

( 3 • le ) 

l'(x0)  = 

3  ( x0 

-  rx)  + 

OJ 

u 

1 

o 

X 

(3. If) 

|l"(xo) 

=  1 

(3.1g) 

if 

(3 . le) 

,  (3. If) 

,  and 

(3.1g) 

are  substituted 

in  (3-ld) 

,  (_3. la)  is 

obtained.  This  was. to  be  expected  because  all  the  exact  developments 
of  qCx)  have  to  be  equivalent. 

The  way  of  developing  q(x)  in  (3. Id)  is  used  in  the  derivation 
of  the  formulas  for  the  application  of  the  proposed  method  to  the 
load-flow  problem,  because  it  yields  the  final  equations  in  less 
steps  than  does  the  Taylor  series  expansion.  This  method  is 
implemented  in  Appendix  1  as  well. 

Equation  (3- la)  can  be  rewritten  in  a  different  way  (3.1h), 
from  which  an  iterative  procedure  for  finding  x%  can  be  derived, 

(x,  -  x0){q'(x0)  +  (x,  -  xQ)}  =  -  q(x0)  ( 3 . lh ) 


and  the  iterative  technique  would  be  (3.2a), 
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xi+l  =  x~  - 


q(xQ) 


o 


q'(x0)+  -  (x.  _  ) 


(3.2a) 


For  the  solution  value  x„ ,  this  gives  (3.2b) 


x*  =  xo  - 


q.U0) 


q'(xo)  +  .  (x,  -  x0) 


(3.2b) 


The  convergence  requirements  for  an  iterative  formula  such  as 
(3*3a) ,  are  the  following: 

Given  the  iterative  sequence, 

xi+l  =  g(xi)  i  =  0,1,2,...  (3.3a) 

whose  initial  estimate  will  be  called  x0 ,  then, 

-  if  l)  g(x)  -satisfies  the  Lipschitz  condition 

|g(x)  -  g(x')|  -  A|x  -  x'|  ,  0  -A<1  (3.3b) 

for  all  the  values  x,  x'  in  the  closed  interval 
I  s  [xQ  -  p,xQ  +  p] ,  where  A  is  the  Lipschitz  cons¬ 
tant;  and 

-  if  2)  the  initial  estimate  xQ  is  such  that 

|xQ  -  g(xQ)|  -  (1  -A  )p  (3-3c) 


then,  a)  all  the  iterates  defined  by  x-{_+p=:g(x^)  lie  within  the 


interval  I, 
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b)  the  iterates  x^  ,  converge  to  some  point 


lim  x-[  =  a 


(3.3d) 


l-yx> 


which  is  a  root,  i.e.:  a  -  g(a);  and 


c)  ot  is  the  only  root  in  the  interval  I 


Furthermore ,  it  can  be  shown  that : 


-if  3)  for  any  x  within  the  interval  I 


g' (x)  |  -  X 


(3.3e) 


and 


-  if  4)  the  initial  estimate  xQ  satisfies  the  condition  2), 


(3-3c)  ,  then, 


Conclusions  a),  b),  and  c)  are  valid. 

Assuming  that  xQ  complies  with  4),  the  convergence  of  the  new  method 
can  be  checked  by  means  of  the  condition  3),  equation  (3-3e). 

During  the  first  iteration,  xQ  is  a  variable,  but  during  the 
other  iterations,  xQ  is  a  constant  parameter.  Therefore,  there  will 
be  a  different  function  g'(xi)  for  the  first  and  the  other  iterations. 
Deriving  (3.2a)  with  respect  to  xQ ,  we  get  g'(x0).  Thus, 


(3.4a) 


and  thus  when  xQ  belongs  to  the  interval  I,  the  initial  value  must 
satisfy, 


S'  (x0)  |  =  x 

<r  (x0)a 


(3.4b) 
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For  the  iterations  beyond  the  first,  xQ  is  regarded  as  a 
constant.  Taking  the  derivation  of  the  right  hand  side  of  (3.2a) 
with  respect  to  xq  gives  (3.4c) 

I  r  \  <l"(xo)  | 

Ig'Upl  =  - 9  X°  ...g- -  <  A  (3.4c) 


{q'(xQ)  +  ^L(Xi  -  x0)}2 


By  subtracting  (3- 2a)  from  (3 -2b),  i.e.:  xq+q  from  x^  an 

expression  for  the  error  at  the  (i+l)  iteration  is  obtained. 


x*  "  xi+l  =  ei+l  = 


q(Xo)  £. 


{q'(x0)  +  (xi  -  xo)}{a'(xo)  + 


dM(x0) 

2 


( x 


(3.5a) 


However,  for  the  first  iteration,  the  expression  is  different. 
Considering  (3.1a),  if  the  expression  is  divided  by  q'(x0),  this 
gives, 


0  = 


lUo) 

q'  (x0) 


+  x 


x 


o 


+ 


<l"(x0)  / 
2qTl^') 


(3-5b) 


but  according  to  (3.2a)  applied  for  xq. 


x 


1 


=  x. 


lUo) 

i'Tx^T 


(3.5c) 


From  (3.5b)  and  (3. 5c),  and  bearing  in  mind  that  £q  -  x*  -  Xq 
and  Zr,  =  x,  -  x„,  the  error  at  the  end  of  the  first  iteration  is 

U  vfc  Q  7 


(3.5d), 
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ei  -  - 


q.n(x0) 
2q' (x0) 


(3.5d) 


A  new  expression  can  be  found  for  Multiplying  the  nume¬ 

rator  and  denominator  of  (3. 5a)  by  £0=(x^  -  x  ),  and  recalling 
(3.1h),  the  expression  (3-5e)  can  be  written, 


ei+l 


q.(x  )qn(x0)  P 

2  o  i 

q(x0){q'(x0)  +  5  (Xi  -  xq ) } 


(3.5e) 


Dividing  the  numerator  and  the  denominator  of  (3-5e)  by 
q(x0)q'(x0),  and  nothing  that  x^_  -  x0  =  £0  -  ,  this  gives, 


-i+1 


1  + 


1  ( X0 ^  F  F  • 

~ - 77 - 7 

2q  1  (x0) 

l" (xo )  ( 

'2qTXx^T  U° 


(3-5f) 


and  taking  the  expression  (3-5d)  into  account, 


'i+1 


£1 

£o 


1  -  fl  +  fi 


'O 


£o  e 


£i_ 

o 


(3.5g) 


For  a  convergent  problem,  when  the  iterates  are  close  to  the 
solution,  £i/e0  will  be  close  to  zero,  and  the  term  £-^£^/eo£0  , 
in  the  denominator  of  (3*5g)j  can  be  dropped.  The  successive  errors 
in  the  neighbourhood  of  the  solution  will  satisfy  the  approximate 
formula  (3-5h) 


s 

e  i+1  “ 


£1 

Eo" 


i 


( 3 . 5h) 


- 
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where  the  subscript  "s"  in  the  error,  indicates  that  these 
errors  are  small  in  comparison  with  the  initial  error,  bacause  the 
iterates  are  close  to  the  solution. 

By  comparing  (3-5d),  (3*5g)5and  (3.5h),  it  can  be  concluded 

that  the  convergence  goes  from  quadratic  at  the  first  iteration  to 
linear  at  the  end  of  the  process. 


3.1.2  THE  NEWTON -RAPHS ON  METHOD  IN  THE  ONE-DIMENSIONAL  CASE.  Con¬ 
sidering  the  development  (3.6), 

q(x#)  =  0  =  q(xi)  +  (x*  -  xi)q’(xi)  +  ^  (x#  -  xi)2q"(xi)  (3-6) 

and  dropping  the  last  right  hand  side  term  of  (3-6),  the  Newton- 
Raphson’s  iterative  procedure  is  obtained, 


xi+l  =  xi  " 


q(y-i) 
q'  (xi) 


(3.7) 


Recalling  the  iterative  formula  (3.3a)  and  the  convergence  con¬ 
dition  (3.3e),  it  can  be  written  for  the  Newton-Raphson  method. 


g' Up) 


q(xj)q,,(xi)  |  < 


l’  Up) 


2 


(3.8) 


Dividing  (3.6)  by  q'(x),  the  error  formula  (3-9)  is  obtained. 


ei+1  =  -  ey 

2q' (xi) 


(3-9) 


which  means  that  the  convergence  is  always  quadratic. 


■ 


. 
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3.1.3  COMPARISON  OF  THE  CONVERGENCE  OF  THE  NEW  METHOD  WITH  A  ONE¬ 
DIMENSIONAL  QUADRATIC  SYSTEM,  WITH  THAT  OF  NEWTON-RAPHSON ’S .  For 
the  Newton-Raphson  method  there  always  exists  an  interval  where 
the  convergence  condition  (3-8)  is  fullfilled,  provided  that  g’(x) 
is  a  continous  function,  "because  g'(x^)  =  0  ,  as  q(x#)  =  0  • 

Since  the  formuals  (3-5d)  and  (3.9)  are  equivalent,  it  can 
"be  concluded  that  the  convergence  of  both  methods  are  exactly 
alike  during  the  first  iteration. 

For  the  convergence  of  any  two  succesive  iterations  of  the 
Newton-Raphson  method,  the  following  condition,  derived  from  (3-9), 
must  be  satisfied  for  the  initial  estimate  xQ,  and  all  the  iterates 


+~i  |  _ 

£ 


q”(xj )  e. |  < 


i  -  l 


i=  0,1,2,... 


-  ■]_  2q'(xi) 

Similarly  for  the  proposed  method,  if  any  ei+i/ei  has  to  be 
equal  or  less  than  one,  from  (3*5g),  it  can  be  derived  that. 


(3.10) 


-  1  * 


^1 


*  1 


(3.11a) 


1  -  fcl 


:1  ei 


eo  £o  £o 


.  The  condition  (3.11a)  must  be  satisfied  to  guarantee  the  con¬ 
vergence  of  the  method. 

The  term  £^/e0  for  every  "i"  ,  can  be  put  in  terms  of  eq/e0  by 
means  of  the  equation  (3-5g)-  The  expression  (3.11a)  with  £i/e0  in 
terms  of  £^/e Q,  becomes  more  and  more  complicated  as  "iM  increases. 
Although  a  general  proof  for  any  "i"  has  not  been  derived,  it  has 
been  found  that  (3.11a)  is  satisfied  for  the  values  of  "i” ,  1,2,3, 
and  U ,  provided  that  ]  e-p /  eQ  J  is  less  than  one. 


■ 
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3.2  DEVELOPMENT  OF  THE  METHOD 

The  proposed  method  can  he  put  in  the  same  classification  as 
the  Newton-Raphson  method;  both  are  direct  methods,  because  the 
solution  of  a  system  of  nonlinear  equations  gives  a  set  of  voltage 
corrections.  As  for  the  Newton-Raphson  method,  the  proposed  method 
is  based  on  the  set  of  nonlinear  equations  (1.5),  thus  the  bus  ad- 
mitance  form  is  used. 

If  the  system  of  nonlinear  equations  (1.5)  is  examined, 

PP  =  ^p^OpPl  +  fqBPl)  +  fp(fqGPl  "  eqBpq) } 

(1.5) 

Qp  =  |(fp(eqGpq  +  fqBpq)  -  ep(fqGpq  -  eqBpq)} 

it  can  be  seen  that  each  equation  consists  of  a  sum  of  terms,  each 
being  a  quadratic  expression  of  the  bus  voltages.  The  general  form 
of  the  terms  is  always  one  of  the  following;  Gijeifj,  Gpjeifj, 

GiJfifjj  Bij  eiej  ,  Bijepfj,  Bijfjfj. 

The  proposed  method  is  based  on  an  iterative  method  of  solution 
of  a  system  of  quadratic  nonlinear  equations  that  is  equivalent  to 
the  system  (1.5).  The  way  to  obtain  the  equivalent  system  and  its 
solution  is  analogous  to  that  of  appendix  1,  where  a  two  dimensional 
example  of  the  procedure  is  shown. 

Let  the  set  e*i,  f*i,  (i  =  1,2,...,N),  be  the  solution  of 
(1.5),  and  let  e°i ,  f°i 5  be  a  set  of  values  which  are  not  the  solu¬ 
tion.  Let  Ae°-j_ ,  Af°-j_  be  the  set  of  increments,  such  that  they 
bring  e°i ,  f°i,  to  the  exact  solution  of  (1.5),  e*£,  f*^. 


, 


(3.13) 


*±  =  e°.  +  Ae° 
i  =  f°i  +  Af° 


substituting  (3.13)  in  (1.5), 


PP  “  q^e°P  +  Ae°p)((e°q  +  Ae°q)Gpq  +  (f°q  +  Af°q)Bpq3  + 

+  ( l°p  +  Af0p)[(f°q  +  Af°q)Gpq  -  (e°q  +  Ae°q)Bpq]} 

(3.14a) 

Qp  =  £{(f°p  +  AfOp)[(eoq  +  Ae°q)Gpq  +  (f°q  +  Af°q)Bpq]- 

-  ( e°p  +  Ae°p)[(f°q  +  Af°q)Gpq  -  (e°q  +  Ae°q)Bpq] } 

calling , 

P°p  -  ^(e°p(e°qGpq  +  P°qBpq)  +  f°p(f°qGpq  -  e°qBpq)} 

Q  p  q^°p(e°q^pq  +  B°qBpq^  “  e°p(BO(lGpq  “  e°qBPQ.)  ^ 

(3.14b) 

and  Pp  =  Pp  -  P°p  ,  AQp  =  Qp  -  Q°p,  then 

Pp  -  P  p  =  E{e°p( Ae°qGpq  +  Af°qBpq)  +  f°p(Af°qGpq  -  Ae°qBpq)}+ 
+  Ae°pZ { ( e°q  +  Aeoq)Gpq  +  (f°q  +  Af°q)Bpq}  + 

+  P0p4{(f°q  +  Af°q)Gpq  -  (e°q  +  Ae°q)Bpq} 

(3.14c) 

QP  -  Q°P  =  |{f°p(Ae°lGPl  +  Af°lBpq)  “  e°p(Af°CLGpq  "  Ae°qBpq) }+ 

+  Ae°p^(e°q  +  AeOq)Bpq  -  (f°q  +  Af°q) Gpq} 

+  Af°pl{(e°q  +  Ae°q)Gpq  +  (f°q  +  AfOq)Bpq} 


+ 
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which  is  equivalent  to 


q^e°q(  e°p^pq  ^°p®pq)  +  Af°q(e°pBpq  +  f’°pGpq)} 

+  Ae°pZ{(e°q  +  Ae°q)Gp^  +  (f°q  +  Af°q)Bpq} 

+  Af°pI{(f°q  +  Af°q)Gpci  -  ( e°q  +  Ae°q)Bpq} 

(3.lUd) 


AQp  =  3{Ae°q( e°pBpq  +  f°pGpq)  +  f°q(f°qBpq  -  e°pGpq)} 

+  Ae°p3 { ( e°q  +  Ae°q)Bpq  -  (f°q  +  Af°q)Gpq) 

+  Af°p2{(e°q  +  Ae°q)Gpa  +  (f°q  +  Af°q)Bpq} 

The  system  (3.1^d)  is  equivalent  to  system  (1.5)*  This  implies 
that  if  the  system  of  nonlinear  equations  (3.1^d)  is  solved,  the 
solution  of  the  load-flow  problem  will  he  known.  The  unknowns  of 
the  system  (3-l^d)  are  the  set  of  increments  Ae°p ,  Af°p. 

It  is  possible  to  arrange  the  equations  (3-l^d)  in  such  a  way 
that  the  subscripts  of  the  left  hand  side  constants  are  arranged  in 
increasing  order;  first  for  APp,  and  afterwards  for  AQp,  as  in 
A?q,  AP2,*-*,  AP^,  AQq,  AQ2,...,  AQj\j.  It  is  also  possible  to 
arrange  the  right  hand  side  terms,  so  that  the  unknowns  have  their 
subscripts  in  increasing  order.  At  first  for  Aep,  and  afterwards 
for  Afp,  as  in  Aeq,  Ae2,»--5  Aejj,  Afq,  Afp,...,  Afjj. 

Assuming  that  the  ordering  of  equations  and  terms  described 
above  has  been  done,  and  writing  (3-l^d)  in  matrix  form,  we  can 
observe  that  the  coefficients  are  all  constant  except  those  in  the 


. 


diagonals  of  the  four  similar  sutmatr ices  into  which  the  matrix 
of  coefficients  can  he  subdivided.  The  coefficients  in  the  four 
diagonals  have  a  constant  part,  which  belong  to  the  sum  which  is 
the  first  term  in  the  right  hand  side  of  the  equations  (3.1^d). 

And  they  have  a  variable  part,  which  depends  on  the  value  of  the 
solution  Ae0^,  Af°-j_.  This  part  is  called  the  variable  part 
because  it  will  change  value  during  the  succesive  iterations,  as 
we  find  better  values  for  the  solution  Ae0^  ,  Af°^. 

The  iterative  method  of  solution  that  can  be  applied  is  the 
same  as  that  in  Appendix  1.  The  nonlinear  system  (3.1^d)  is 
linearized  by  giving  an  approximate  value  to  the  unknown  part  of 
the  variable  part  of  the  diagonal  coefficients.  For  instance,  in 

the  diagonal  coefficient  for  Ae°p,  the  variable  part  is, 

E{(e°q  +  Ae°q)Gpq  +  (f°q  +  Af°q)Bpq}  (3.15a) 

since  the  solution  of  (3.1^d)  is  part  of  the  variable  part  of  the 
coefficient,  and  the  solution  is  not  known,  a  value  for  Ae°-j_ ,  and 
Af0^,  will  have  to  be  assumed.  In  the  first  iteration,  these  value 
will  be  taken  as  zero.  In  other  words,  Ae°^ ,  and  Af°^  will  be 
dropped  before  e0^,  and  f°j_.  The  variable  part  of  the  diagonal 
coefficient  of  Ae°p  will  be  then, 

+  fvw  (3-i5b) 

which  is  constant.  Doing  the  same  with  all  the  diagonal  coef¬ 
ficients,  a  linearized  system  results.  Its  solution  will  be  the 
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set  of  values  Ae^-j_ ,  Af^j_.  These  values  can  he  substituted  in 
(3.15a).  The  variable  part  of  the  diagonal  coefficient  of  Ae°p, 
is  then, 


E{(e°n  +  Ae1q)Gpq  +  (f°n  +  AfijB^} 


G'  P<1J 


(3.15c) 


which  is  again  a  constant.  Doing  the  same  with  all  the  diagonal 
elements  ,  a  new  linearized  system  of  equations  will  be  obtained. 
Once  solved,  it  will  give  a  new  set  of  solutions  Ae^,  Af^q. 

It  can  be  expected  that  Ae^j_ ,  Af^-j.  is  a  better  approximation 
to  Ae°j_ ,  Af  °£ ,  than  Ae^ ,  Af^,  because  (3.15c)  is  closer  to 
(3.15a)  than  (3.15b)  is.  The  same  process  can  be  applied  again, 
i.e.:  substitute  Ae  ,  Af^  in  (3.15a)  and  obtain  a  new  and 
better  set  Ae^,  Af^. 

When,  at  the  first  iteration,  Ae°p,  Af0p  are  dropped  before 
e°p,  f °p ,  equation  .(3. lUd)  becomes, 


APp  -  Z{Aeq(e°pGpq  -  f°pBpq)  +  Afq  (e°pBpq  +  f°pGpq)}  + 


+  ^epq^e°Q.^P(l  +  ^0q.®P(l^  + 


+  AfpX{f°<1Gpq  -  e0qBpq} 


(3.16a) 


AQp  “  £(Aeq(e°pBpq  +  ^°P^pq)  +  Afq(f°pBpq  e°p^pq)^  + 


+  Aer)S{e°nB-Da  -  f°qGpq}  +  AfpE{e°qGpq  +  f°qBpq} 


'PqLC  q/GPq. 


(3.16a)  would  give  the  set  Ae1^ ,  Af1^  as  a  solution.  In  general  at 
the  (k+1)  iteration  the  linearized  system  (3.l6b)  would  be  solved 
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App  =  Z{Aeq(e°pGpq  -  f°pBpq)  +  Afq(e°pBpq  +  f°pGpq)>  + 

+  AepZ{ ( e°q  +  Aekq)Gpq  +  (f°q  +  Afkq)Bpq}  + 

+  ^Pq^^  q.  +  Afkq)Gpq  -  (e°q  +  Ae  q)Bpq} 

(3.16b) 

AQp  =  Z{Aeq(e°pBpq  +  f°pGpq)  +  Afq(f°pBpq  -  e°pGpq)}  + 

+  AepZ{(e°q  +  Aekq)Bpq  -  (f°q  +  Afkq)Gpq}  + 

+  AfpZ{(e°q  +  Aekq)Gpq  +  (f°q  +  Afkq)Bpq} 

which  would  yield  the  set  of  values  Aek+'*'i,  Afk+^£,  as  solution. 

3.3  THE  BASES  FOR  CONVERGENCE 

Observing  (2.7),  and  (3.16a)  in  the  former  section,  it  can  be 
seen  that  they  are  the  same.  (2.7)  is  the  general  linearized  system 

which  is  solved  at  each  iteration  of  the  Newton-Raphson  method, 
and  (3.l6a)  is  the  linearized  system  which  is'  solved  at  the  first 
iteration  of  the  proposed  method.  From  this  it  can  be  concluded 
that  the  first  iteration  of  the  proposed  method  is  exactly  the 
same  as  that  of  the  Newton-Raphson  method,  provided  that  the  same 
starting  value  for  the  voltages. is  taken. 

Therefore,  as  far  as  the  first  iteration  is  concerned,  the 
same  reasoning  can  be  applied  as  for  the  numerical'  properties  of 
the  Newton-Raphson  method,  this  means  that  if  the  set  of  conditions 
specified  in  section  2.2  are  satisfied  for  the  initial  guess  of 
the  voltages,  then  the  new  voltages,  which  result  from  correcting 
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the  initial  guess  with  the  increments  that  result  from  the  first 
iteration  are  always  closer  to  the  solution  than  the  voltages  of  the 
first  guess  are.  Moreover,  in  the  case  of  the  proposed  method,  it 
means  as  well  that  all  Ae^i,  Af ^ ,  point  towards  the  solution. 

As  it  has  been  shown  in  section  3.2,  the  variable  parts  of  the 
diagonal  coefficients,  once  the  values  of  Ae^-j_ ,  Af^,  have  been 
substituted,  are  closer  to  the  true  values  of  the  diagonal  coef¬ 
ficients  than  that  of  the  first  iteration.  The  non-diagonal  coef¬ 
ficients  in  the  proposed  method  are  unchanged  from  the  first  until 
the  last  iteration. 

As  for  the  one  dimensional  case  in  section  3.1,  “the  results 
in  the  second  iteration  Ae^-j_ ,  Af  ,  will  be  closer  to  the  solution 
Ae°i  j  Af°i5  than  the  set  Ae-^,  Af^-j_ .  This  is  so  because  the 
linearized  system  (3.16b)  for  iteration  2  is  closer  to  the  exact 
system  (3-l^d)  than  the  linearized  system  (3.16a)  for  the  first 
iteration  is. 

In  other  words,  since  the  non-diagonal  coefficients  are  the 
same  for  the  exact  system  and  for  the  approximate  system,  it  can 

be  expected  that  if  the  value  of  the  diagonal  coefficients  tend 
to  their  true  value,  the  solutions  at  each  iteration  have  to  tend 
towards  the  true  value  of  the  solution  of  the  linear  system  (3.l4d). 

Calling  AS  the  vector  of  differences  of  powers,  T  the  voltage 
vector,  L  the  matrix  of  coefficients  expressing  the  linear  part  of 
the  system  of  equations  (3.1^b),  and  D  the  matrix  of  nonlinear  dia¬ 
gonal  elements,  (3.l4b)  can  be  written  in  the  simplified  way  (3.17a), 


i  j 
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AS  =  [L  +  D(T) ] (T  -  T0) 


(3.17a) 


from  which,  the  iterative  expression  (3.17b)  is  obtained 


Ti+i  =  T0  +  [L  +  D(Ti)]  XAS 


(3.17b) 


Applying  the  condition  (3*3b),  a  sufficient  convergence  con¬ 
dition  (3.18),  can  be  written 


II  jfTtL  +  D(Ti)  ]-1AS  ||  <  1 


(3.18) 


For  the  first  iteration,  the  same  sufficient  convergence  con¬ 
dition  as  for  the  Newton-Raphson  method,  can  be  applied.  Therefore, 


the  requirements  described  in  section  (2.2)  apply  for  the  first 


iteration  of  the  proposed  method. 

The  experimental  results  with  the  IEEE  l4,  30,  and  57  bus 
test  systems  confirm  these  assertions. 

3.4  COMPARISON  BETWEEN  THE  PROPOSED  METHOD  AND  THE  NEWTON-RAPHSON 


METHOD 


Although  the  differences  between  the  Newton-Raphson  method 
and  the  proposed  method  have  been  indicated  already,  more  insight 
could  be  obtained  if  some  of  the  points  are  stressed,  and  some 
of  the  parts  of  both  methods  are  put  side  by  side. 

For  the  sake  of  clarity  while  the  comparison  of  the  two 
methods  is  carried  out ,  the  increment  with  a  point  inside  "A"  will 
be  used  to  refer  to  the  increments  that  appear  in  the  proposed 
method . 

Let  CIRp  and  Clip  be  the  real  and  imaginary  part  of  the  current 


' 

flowing  into  the  system  at  bus  "p".  It  is  known  that  when  the  bus 
voltages  are  e*-j_,  f#^»  the  real  and  imaginary  component  of  the 
current  at  bus  "p”  are, 


CIR*P  =  !<e*qGpq  +  f*qBpq> 

CII*p  =  |(f*qBpq  -  e*qBpq)  (3.19) 

oberving  (2.7)  it  can  be  seen  that  (3.19)  could  be  substituted  in 
there.  A  new  expression  of  (2.7)  could  be  written, 

^Pp  ~  ^^eq(ep9pq  -  fpBpq)  +  Afp(epBpq  +  fpGpq)}  + 

+  AepCIRp  +  Af  Clip 

AQp  =  ^^eq(ep®pq  +  Ip^'pq)  +  ^Pq(^p®pq  “  ePGPl^  + 

+  AfpCIRp  -  AepCIIp  (3o20) 


As  in  (3.19)5  e°-j_ ,  f°-j_,  being  the  initial  guess  voltages,  it 
is  possible  to  define, 


CIR°P  =  ^(e°lGPl  +  f°qBpq) 


CII°  =  ^(f°nB 


q^pq 


-  e°nB  ) 


q  pq; 


(3.21a) 


Let  DCR°p  and  DCI°p  be  the  real  and  imaginary  part  of  the 
difference  between  the  current  at  bus  "p”  when  the  voltages  at  the 
nodes  are  the  solution  voltages,  and  the  current  at  bus  "p"  when 
the  voltages  are  the  initial  guess, 
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DCROp  =  |(AeO(1Gpq  +  AfOqBpq) 

CH°P  =  £(Af°qBpq  "  Ae<VW  (3.21b) 

The  expressions  (3.21a)  and  (3.21b)  can  be  substituted  in 

(3.l4d)  to  get 

=  3{Ae°q( e°pGpq  -  f°pBpq)  +  Af°q(e°pBpq  +  f°pGpq)}  + 

+  Ae°p(CIR°p  +  DCR°p )  +  AfOp(CIlop  +  DCI°p) 

AQp  =  Z{Ae°q(e°pBpq  +  f°pGpq)  +  Af°q(f°pBpq  -  e°pGpq)}  + 

+  Af°p( CIR°p  +  DCR°p)  -  Ae°p(CII°p  +  DCI°p)  (3.22) 

(3.22)  is  an  equivalent  system  to  (3.1^d).  It  can  be  noted  that  in 
the  proposed  method  the  recalculation  of  the  new  diagonal  coefficients 
at  each  iteration  is  equivalent  to  finding  a  better  set  of  incre¬ 
ments  of  currents  DCR-j_,  DCI^. 

Obviously  when  at  the  first  iteration  Ae°p>  Af°p>  are  dropped 
before  e°p,  f°pJ  it  is  the  same  as  dropping  DCR°p ,  DCI°p,  before 
CIR°p ,  CII°p. 

3 - U . 1  DIFFERENCE  BETWEEN  THE  CONCEPT  OF  INCREMENT  OF  VOLTAGE  IN 
THE  NEWTON-RAPHSON  METHOD  AND  IN  THE  PROPOSED  METHOD.  The  difference 
between  the  Newton-Raphson  method  and  the  proposed  method  lies 
basically  in  their  meaning  of  the  increment  of  voltage. 

In  the  Newton-Raphson  method  the  increments  of  voltage  Ae-j_ , 


' 
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are  the  difference  "between  the  values  of  the  voltages  at  two  succe- 
sive  iterations.  In  the  proposed  method  the  increments  of  voltage 
Aej.,  Afi,  are  the  difference  between  the  voltage  calculated  at  the 
end  of  the  iteration,  and  the  first  assumed  value  of  the  voltages. 

Hence,  whereas  the  Newton-Raphson ’ s  increments  of  voltage  be¬ 
come  smaller  and  smaller  at  each  iteration,  the  proposed  method  has 
increments  Aei 5  that  tend  towards  Ae°£,  Af0^,  which  are  different 

from  zero. 

Comparing  (3. 20)  and  (3.22),  it  can  be  seen  that  the  variable 
terms  in  the  variable  part  of  the  diagonal  coefficients  in  the 
Newton-Raphson  method  are  the  bus  currents  calculated  at  the  itera¬ 
ted  voltages,  and  are  the  same  for  the  proposed  method,  because. 


CIR°p  +  DCRp  =  CIRp 
CII°p  +  DCIp  =  Clip 


(3.23) 


It  must  be  kept  in  mind  though  that  the  constant  part  of  the  diago¬ 
nal  coefficients  is  variable  at  each  iteration  in  the  Newton-Raphson 
method,  but  is  constant  in  the  proposed  method. 

One  more  fact  can  be  observed,  APp,  AQp,  are  the  differences 
between  the  scheduled  powers  and  the  powers  computed  with  the  ite¬ 
rated  voltages,  since  the  iterated  voltages  change  at  each  iteration 
APp,  AQp  will  be  different  too,  and  as  the  iterated  voltages  tend 
towards  the  true  voltages,  the  computed  powers  will  tend  to  the 
scheduled  powers,  so  APp,  AQp,  tend  to  zero. 


' 


In  the  proposed  method  APp ,  AQp  are  constant  because  they  are 
the  difference  between  the  scheduled  powers  and  the  powers  calcu¬ 
lated  with  the  first  assumed  value  of  the  voltages. 

3.4.2  CALCULATION  OF  THE  ELEMENTS  OF  THE  JACOBIAN.  From  an  exami¬ 
nation  of  (2.7)  and  (3.l4d),  it  can  be  noted  that  whereas  in  (3*l4d) 
the  non-diagonal  coefficients  are  constant,  because  they  depend  on 
the  value  of  the  first  estimate  of  the  voltages,  the  non-diagonal  coef¬ 
ficients  of  the  Jacobian,  when  the  Newton-Raphson  method  is  used, 
depend  on  the  value  of  the  iterated  voltages,  so  they  have  to  be  re¬ 
calculated  at  each  iteration. 

This  is  one  of  the  most  important  assets  of  the  proposed  method. 

As  far  as  computer  time  is  concerned,  the  recalculation  of  the 
elements  of  the  linearized  system  is  the  largest  item.  This  is  so 
because  of  the  large  size  of  the  matrix  of  the  linearized  system, 

2(NB  -  S)x2(NB  -  S),  (NB  number  of  buses,  S  number  os  slack  buses), 
and  because  of  the  logical  instructions  involved  in  the  handling 
of  the  elements  of  the  bus  admittance  matrix. 

If  with  the  proposed  method,  advantage  is  to  be  taken  of  the 
fact  that  the  non-diagonal  coefficients  are  the  same  throughout  all 
the  iterations,  a  non-destructive  technique  for  solving  the  linear 
system  (3*l6b),  must  be  used.  If  a  destructive  technique  is  imple¬ 
mented,  extra-storage  must  be  provided  to  save  the  constant  elements 


of  the  matrix. 


. 

* 
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3*5  OUTSTANDING  FEATURES  OF  THE  PROPOSED  METHOD 

As  a  brief  outline  of  the  characteristics  of  the  proposed  method, 
it  could  be  said  that: 

-  The  proposed  method  is  a  direct  method  operating  with  the  bus 
admittance  form. 

-  The  proposed  method  implements  a  new  way  to  solve  the  system 
of  nonlinear  equations  (1.5),  which  is  suitable  for  a  system 
of  quadratic  nonlinear  equations. 

-  In  the  proposed  method  the  only  coefficients  that  change  with 
the  iterated  voltage,  are  the  diagonal  ones,  so  they  have  to 
be  recalculated  at  each  iteration. 

-  The  convergence  conditions  of  the  proposed  method  are  diferent 
from  those  of  the  Newton-Raphson  method.  The  conditions  for 
the  one-dimensional  case  are  in  section  3*1.3  * 

-  During  the  first  iteration  the  rate  of  convergence  of  the 
proposed  method  is  quadratic.  The  rate  of  convergence  is  less 
than  quadratic  and  decrease  as  the  number  of  iterations  increase. 
The  rate  of  convergence  becomes  linear  in  the  neighbourhood  of 
the  solution. 

-  It  is  very  easy  to  switch  from  the  proposed  method  to  the 
Newton-Raphson  method  and  back,  as  far  as  computer  programming 
is  concerned.  A  flow  chart  of  the  proposed  method  can  be 


found  in  Appendix  3* 


IV. 


THE  PROPOSED  METHOD  USING  TRIANGULATION 


WITH  PIVOTING  FOR  SOLVING  THE  LINEARIZED 

SYSTEM  AT  EACH  ITERATION.  RESULTS 
AND  CRITICISM 

4.1  NUMERICAL  RESULTS  OF  THE  PROPOSED  METHOD 

In  this  chapter  the  performance  of  the  proposed  method  using  triangul¬ 
ation  with  pivoting  is  examined.  The  experimentation  done  has  shown  that 
other  ways  than  the  exact  solution  by  triangulation  with  pivoting  are  more 
suitable  for  the  solution  of  the  linearized  system  (3.lUd)  at  each  iteration 
when  using  the  proposed  method. 

Nevertheless,  since  triangulation  with  pivoting  is  the  usual  technique 
for  the  Newton-Raphson  method,  it  is  worth  knowing  what  the  proposed  method 
does  with  it,  and  analyzing  the  results  obtained. 

4.1.1  EXPLANATIONS  ON  THE  COMPUTER  PROGRAM  USED.  The  goal  of  this  chapter 
is  the  comparison  of  the  methods  using  the  technique,  of  solution  of  a 
linear  system  of  equations,  usually  implemented  by  the  Newton-Raphson 
method,  namely,  triangulation  with  pivoting. 

In  order  to  take  advantage  of  the  constancy  of  the  non-diagonal  coef¬ 
ficients  of  the  linearized  system,  double  computer  storage  must  be  used 
for  the  coefficients  because  triangulation  with  pivoting  destroys  the 
linear  system  when  it  finds  its  solution. 

In  the  computer  program  used,  the  whole  bus  admittance  matrix  is 

used,  and  no  logic  programming  complications  exist  for  handling  the 
elements  of  this  matrix.  This  is  pointed  out  because  in  other  sections 


■ 
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the  load-flow  program  used  implements  sets  of  logic  instructions 
aimed  at  avoiding  the  storage  of  the  zero  elements  of  the  bus  ad¬ 
mittance  matrix. 

Whenever  a  program  which  measures  load-flow  solution  times 'has 
been  run,  it  has  been  run  along  with  a  reference  case  so  that  the 
times  can  always  be  compared.  It  is  known  that  a  computer  can  spend 
different  amounts  of  time  on  the  same  problem,  due  to  the  different 
modes  of  shared  time. 

There  is,  in  the  load-flow  problem  shown  in  Appendix  3,  a  trial 
and  error  technique  that  takes  into  account  the  constraints  in  the 
available  reactive  powers  at  the  voltage  controlled  buses.  However, 
for  measuring  the  convergence  and  speed  of  solution  of  the  load-flow 
method,  the  first  run  of  the  trial  and  error  technique  is  enough. 

This  is  what  has  been  made  throughout  this  thesis. 

4.1.2  TEST  SYSTEM  USED.  The  IEEE  14  and  30  bus  test  systems  were 
used  for  making  a  comparison  between  the  Newton-Raphson  method  and  the 
proposed  method.  The  results  have  been  different  as  far  as  time  is 
concerned  between  the  two  methods,  but  similar  for  each  test  system. 

For  the  sake  of  brevity,  only  the  results  obtained  with  the  30 
bus  test  system  will  be  shown. 

4.1.3  NUMERICAL  RESULTS.  The  tolerance  used  in  the  comparison  ranges 
from  0.0001  to  0.0002  units  of  power  or  voltage,  and  will  be  indicated 
in  every  case.  The  precision  will  be  called  EPS. 

As  can  be  seen  in  Appendix  3  there  are  two  ways  to  get  out  of  the 
iterative  procedure  of  the  proposed  method.  One  by  meeting  the  precision 


' 
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in  power  at  the  buses,  and  the  other  when  the  changes  of  voltage  at 
the  iteration  are  less  than  the  required  precision,  EPS. 

It  has  been  observed  that,  in  most  of  the  cases,  when  the  change 
in  voltage  is  less  than  EPS,  the  mismatch  in  power  at  the  buses  with 
respect  to  the  scheduled  powers  is  not  necessarily  less  than  EPS.  An 
error  can  thus  be  found  at  the  end  of  the  process.  It  could  be  elimi¬ 
nated  by  requiring  a  higher  precision  in  voltage  change  to  be  met, 
but  there  are  some  reasons  for  not  doing  so.  A  further  discussion  about 
these  points  can  be  found  in  the  section  where  there  are  the  comments 
on  the  advantages  and  disadvantages  of  the  proposed  method.  In  this 
section  it  will  be  seen  as  well  that  the  rate  of  convergence  decreases 
as  the  solution  is  approached. 

The  performance  of  the  Newton-Raphson  method  and  the  proposed 
method  are  shown  in  Table  (4.1). 


TABLE  4.1  PERFORMANCE  OF  THE  NEWTON-RAPHSON  METHOD  AND  THE  PROPOSED 
METHOD  USING  TRIANGULATION  WITH  PIVOTING;  ON  THE  IEEE  30 
BUS  TEST  SYSTEM. 


Newton-Raphson 


EPS  0.0001 

Complete  iterations  3 

Iterations  just  started  1 

Total  number  of  iterations  4 

Time  for  the  1st  iteration  1.639  sec. 


Time  for  iterations  other  than  the  1st  1.689  sec. 


Proposed  method 
0.0002 
5 
0 
5 

1.684  sec. 
1.476  sec. 


Total  time  to  meet  EPS 


5.127  sec. 


7.589  sec. 


r 
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The  same  starting  values  for  the  voltages  have  been  used  with  each 
method.  These  values  are,  e0^  =  1.0  and  f°^  =  -  0.3.  The  slack  bus 
has  the  values  es  =  1.06  and  fs  =  0.0  . 

The  counting  of  the  time  includes  only  the  time  spent  in  the 
iterations.  The  writing  of  the  scheduled  powers  and  the  initial  volt¬ 
ages,  the  calculation  of  the  line  flows,  and  the  calculations  of  the 
reactive  powers  at  the  voltage  controlled  buses  at  the  end  of  the 
problem,  are  not  included. 

The  proposed  method  using  triangulation  with  pivoting,  conver¬ 
ges,  but  not  as  fast  nor  precisely  as  the  Newton-Raphson  method. 
Comments  on  the  .precision  can  be  found  in  section  4.4  . 

4.2  THE  EFFECT  OF  ACCELERATION  FACTORS 

Acceleration  factors  have  been  tried  with  the  proposed  method 
using  triangulation  with  pivoting  for  solving  the  linearized  system. 

The  explanation  of  the  convergence  acceleration  of  a  load-flow  method 
by  means  of  acceleration  factors,  can  be  found  in  section  8.4  of  (3). 

The  way  in  which  they  have  been  used  is  by  multiplying  the  bus 
voltage  changes  at  each  iteration  by  a  factor.  A  pair  of  different 
factors  were  used,  one  for  the  real  component  of  the  voltage,  and 
another  for  the  reactive.  The  value  of  the  factors  were  near  unity. 

Acceleration  factors  were  tried  in  two  ways:  since  the  first 
iteration,  and  from  the  second  iteration  on. 

The  improvement  observed  with  the  best  pair  of  acceleration  fac¬ 
tors,  was  not  enough  to  get  to  the  solution  with  less  iterations  than 
without  acceleration  factors. 


' 
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4.3  COMBINATION  OF  THE  NEWTON- RAPHS ON  METHOD  AND  THE  PROPOSED  METHOD 
USING  TRIANGULATION  WITH  PIVOTING 

It  can  be  seen  in  appendix  3  that  it  is  very  easy  to  switch  from 
the  Newton-Raphson  method  to  the  proposed  method  and  back.  Taking  now 
advantage  of  this  fact,  the  30  bus  load-flow  problem  can  be  solved  by 
an  alternation  of  the  Newton-Raphson  method  and  the  proposed  method. 

The  method  is  convergent,  faster  than  the  proposed  method  but  not  as 
fast  as  the  plain  Newton-Raphson  method.  The  results  are: 

EPS  =  0.00015 

number  of  iterations  =  4  (4th  iteration  complete) 

time  per  iteration:  for  iterations  1  and  3,  (Newton-Raphson)  =  1.67  sec. 

for  iterations  2  and  4,  (proposed  method)=  1.57  sec. 
total  time  to  meet  EPS  =  6.485  seconds 

Other  combinations  of  the  Newton-Raphson  method  and  the  proposed 
method  were  tried.  In  one  of  them,  a  switch  to  the  Newton-Raphson 
method  is  done  when  the  value  of  the  changes  in  voltage  at  the  ite¬ 
rations  are  less  than  a  specified  value.  The  load-flow  problem  has 
been  solved  with  0.03  as  the  switch  value  for  the  changes  of  voltage. 

The  results  obtained  have  been: 

EPS  =  0.0001 

number  of  iterations  =  4  (4th  iteration  complete) 

times  per  iterations:  for  iters.  1,  3  and  4,  (Newton-Raphson)  =  1.64  sec. 

for  iteration  2,  (proposed  method)  =  1.48  seconds 


total  time  to  meet  EPS  =  6.404  seconds 
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It  must  be  pointed  out  that  the  same  precision  as  for  the  Newton- 
Raphson  has  been  used. 

The  load-flow  problem  has  finally  been  solved  implementing  the 
proposed  method  using  triangulation  with  pivoting  but  using  the  Newton- 
Raphson  method  when  any  of  the  voltage  changes  at  the  former  iteration 
is  greater  than  0.1  .  The  results  are: 

EPS  =  0.0001 

number  of  iterations  =  4  (4th  iteration  complete) 

time  per  iteration:  for  iterations  1  and  2  (Newton-Raphson)  =  1.63  seconds 

for  iterations  3  and  4  (proposed  method)^  1.52  seconds 

total  time  to  meet  EPS  =  6.132  seconds 

4.4  THE  WEAK  POINTS  IN  THE  PROPOSED  METHOD 

By  analyzing  the  weak  points  of  the  proposed  method  using  trian¬ 
gulation,  and  those  of  Newton-Raphson ’ s  in  section  4.5,  the  changes  to 
apply  to  the  proposed  method  can  be  inferred. 

-  The  rath  of  convergence  decreases  as  the  solution  is  approached. 

This  means  that  unless  the  recalculation  of  the  terms  of  the  linearized 

system  takes  a  great  amount  of  time  at  each  iteration,  the  extra  ite¬ 
rations  needed  will  offset  the  savings  in  time  per  iteration. 

Table  2.1  shows  that,  in  the  Newton-Raphson  method,  the  time  for 
the  calculation  of  the  jacobian  ranges  from  47.6 %  for  the  14  bus  case 
to  58.7%  of  the  iteration  time  for  the  57  bus  case.  With  the  30  bus 
case,  the  calculation  of  the  new  jacobian  takes  54.5%  of  the  iteration 
time.  The  reason  why  in  table  4.1  only  a  12.3%  reduction  of  time  per 
iteration  is  observed,  is  that  the  results  of  table  2.1  correspond 
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to  a  program  which  only  stores  in  the  computer  memory  the  nonzero 

elements  of  ,  and  therefore  has  to  have  logical  instructions  to 

work  with  the  non-zero  elements  of  ,  which  are  time  consuming. 

The  terms  of  the  jacobian  are  a  function  of  G  and  B  .  (Y 

pq  pq  pq 

-  JBp^) .  On  the  other  hand  the  results  in  table  4.1  come  from  a 
program  in  which  the  whole  matrix  Y^us  :''s  store<^>  so  no  logical 
instructions  are  needed.  In  section  8.2  the  reasons  for  storing  only 
the  nonzero  elements  of  Y^g  are  discussed. 

-  Another  major  inconvenience  of  the  proposed  method  using 
triangulation  with  pivoting  is  the  extra  storage  required  for  keeping 
the  coefficients  of  the  linearized  system.  The  solution  of  a  system  of 
linear  equations  by  triangulation  with  pivoting  destroys  the  matrix  of 
coefficients  and  the  vector  of  constants.  Therefore,  for  not  having 

to  recalculate  the  whole  set  of  coefficients  at  each  iteration,  double 
storage  must  be  provided.  The  coefficients  of  the  linearized  system 
mean  £ 2 (NB  -  S)]  x  [ 2 (NB  -  S)]  words,  where  NB  is  the  number  of  buses 
and  S  is  the  number  of  slack  buses. 

The  load-flow  problem  with  the  proposed  method  using  triangul¬ 
ation  with  pivoting,  and  storing  only  the  non-zero  elements  of  Y^ug 
has  not  been  tried  because  even  though  the  proposed  method  could  be 
as  good  or  even  better  than  Newton-Raphson  method  as  far  as  time  is 
concerned,  the  extra  storage  for  the  coefficients  of  the  linearized 
system  offsets  any  time  advantage. 

-  Another  of  the  inconveniences  of  the  proposed  method  is  that 
an  iteration  can  be  reached  where  the  voltage  changes  are  less  than 
the  precision  EPS  but  the  bus  powers  are  not  yet  within  EPS 
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neighbourhood  of  the  scheduled  powers. 

Succesive  iterations  of  the  method  lead  the  powers  towards 
the  scheduled  powers.  This  would  mean  a  waste  of  time,  because  a 
very  little  change  in  voltage  would  take  some  whole  iterations  and 
there  is  a  risk  of  getting  entangled  due  to  the  roundoff  error. 

A  measurement  of  the  mismatch  of  power  can  be  found  in  table 
6.10,  where  the  missmatch  is  compared  with  the  changes  in  voltage 

in; the  iterations,  and  an  explanation  is  given. 

A.  5  THE  WEAK  POINTS  IN  THE  NEWTON-RAPIISON  METHOD 

The  disadvantages  that  will  be  described  below  are  also  dis¬ 
advantages  of  the  proposed  method  during  the  first  iteration,  due 
to  the  fact  that  the  first  iteration  of  the  proposed  method  is  as 
Newton-Raphson’ s . 

The  rate  of  convergence  of  the  Newton-Raphson  method  is 
quadratic.  But  it  is  not  worth  spending  so  much  time  on  the  first 
iteration  for  solving  a  linearized  system,  because  beforehand,  we 
know  that  it  is  going  to  yield  only  an  approximate  answer. 

The  time  spent  on  the  solution  of  the  linear  system  ranges  from 
47.3%  of  the  iteration  time  in  the  IEEE  14  bus  test  system,  to  41.2% 
in  the  57  bus  test  system. 

The  Newton-Raphson  method  has  to  recalculate  the  elements  of 
the  jacobian  at  each  iteration,  which  takes  a  long  time. 

4.6  POSSIBLE  IMPROVEMENTS 

The  results  of  the  thesis  show  that  improvements  are  indeed 
possible.  The  characteristics  that  an  alternative  method  should  have 
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to  give  an  improvement  in  time  when  used  with  the  proposed  method  are 
explained  below. 

4.6.1  THE  EFFICIENCY  OF  THE  STEPS  OF  A  METHOD  WITH  RESPECT  TO  TIME. 
Examining  the  proposed  method  it  could  be  seen  that  if  a  non-destruc¬ 
tive  system  of  solution  of  a  linear  system  were  used,  advantage  could 
be  taken  of  not  having  to  recalculate  the  non-diagonal  coefficients 
of  the  linearized  system  at  each  iteration.  Besides  a  double  storage 
for  the  coefficients  would  not  have  to  be  provided. 

The  triangulation  with  pivoting  method  for  the  solution  of  a 
system  of  linear  equations  is  an  exact  method  of  solution.  It  is 
impossible  to  obtain  with  this  method  an  approximate  solution  using 
less  time,  because  once  the  method  is  started,  it  must  be  carried  out 
until  the  end.  Since  at  the  first  iteration  some  of  the  coefficients 
are  not  close  to  the  solution,  an  exact  solution  of  the  linearized 
equation  is  a  waste  of  time.  An  approximate  less  time  consuming 
solution  would  be  more  suitable. 

4.6.2  ADVANTAGE  OF  THE  PROPOSED  METHOD  OVER  THE  NEWTON-RAPHSON  METHOD. 
If  an  approximate  method  of  solution  of  the  linearized  system  is 
implemented  during  the  early  stages  of  the  load-flow,  more  iterations 
will  have  to  be  used  before  arriving  to  the  solution. 

If  the  Newton-Raphson  method  were  employed  along  with  an  approxi¬ 
mate  solution  of  the  linearized  system  during  the  early 
stages,  no  time  improvement  would  be  obtained  because,  at  each 
iteration,  it  would  be  necessary  to  recompute  all  the  terms  of  the 


jacobian,  which  alone  takes  the  58.7%  of  the  iteration  time  in 
the  Newton-Raphson  solution,  using  triangulation  with  pivoting, 
when  applied  to  the  IEEE  57  bus  test  system.  Using  the  proposed 
method  this  inconvenience  would  not  exist  because,  at  each  iteration, 
only  the  diagonal  coefficients  would  have  to  be  recalculated.  This 
time  saving  due  to  calculation  of  the  diagonal  elements  only  can 
be  expected  to  increase  with  the  size  of  the  system. 

A  time  improvement  can  thus  be  expected  with  the  proposed  method 
using  an  approximate  technique  for  the  solution  of  the  linearized 
system. 
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V. 


THE  THEORETICAL  BASIS  FOR  THE  USE  OF  THE 


GAUSS-SEIDEL  ITERATIVE  TECHNIQUE  FOR 

SOLVING  THE  LINEARIZED  SYSTEM 
IN  THE  PROPOSED  METHOD 

5.1  FEASIBILITY  OF  THE  SOLUTION 

In  this  chapter  it  is  shown  how  the  Gauss-Seidel  iterative 
technique  for  solving  a  system  of  linear  equations  complies  with 
the  requirements  pointed  out  in  section  4.6,  and  how  the  propos¬ 
ed  method  coupled  with  the  Gauss-Seidel  technique  can  yield 
results  comparable  with  those  obtained  from  the  Newton-Raphson 
method  using  triangulation  with  pivoting. 

Sections  7.1  and  7.2  of  chapter  7,  illustrate  how  acceleration 
factors  successfully  increase  the  speed  of  the  Gauss-Seidel  technique 

5.1.1  REQUIREMENTS  FOR  CONVERGENCE.  A  number  of  sufficient 
conditions  can  be  found  in  the  literature  of  numerical  analysis, 
for  a  linear  system  to  have  a  convergent  iterative  solution  when 
using  the  Gauss-Seidel  technique.  An  explanation  of  Gauss-Seidel 
or  successive  iteration  technique  can  be  found  in  section  4.2 

of  (11) ,  and  some  tests  of  convergence  are  also  presented  in 
section  4.4  of  (11). 

The  convergence  conditions  are  generally  based  on  the 
fulfilment  of  some  numerical  conditions  by  the  elements  of  the 
matrix  of  coefficients.  Unfortunately  a  general  theoretical 
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proof  of  convergence  would  be  very  long  and  cumbersome,  and  a  numerical 
proof  for  a  specified  system,  is  useless,  because  it  does  not  guarantee 
the  convergence  for  any  other  case. 

In  section  5.9  of  Conte's  book  on  numerical  methods  (12), 
there  is  a  sufficient  but  not  necessary  convergence  condition  for  the 
Gauss-Seidel  technique.  The  use  of  this  condition  does  not  yield  a 
general  convergence  proof,  but  applying  the  particular  condition  to  the 
proposed  method,  provides  us  with  some  insight  into  the  general  conver¬ 
gence  properties  of  the  linearized  system. 

Let  us  denote  Ipq  as  a  general  term  for  the  matrix  of  linearized 
coefficients,  which  will  be  called  L  henceforth.  Then  a  sufficient, 
but  not  necessary  convergence  condition  for  the  Gauss-Seidel  technique 
is  that,  for  any  row  "p", 


1pq|  < ] 1pq| ’  q  ^  p 


(5.1) 


5.1.2  CHARACTERISTICS  OF  THE  MATRIX  OF  COEFFICIENTS  OF  THE  LINEARIZED 
SYSTEM.  Equation  (3.14d)  can  be  written  in  the  simplified  way. 


_AP_ 

_Lf_  ! 

L11  ~ 

X 

4?_ 

;  where 

L1  iL11 

AQ 

Lin  i 

LIV 

Af 

1 

LlII|  LIV 

=  L 


(5.2a) 


The  diagonal  elements  of  each  submatrix  L^,  L^,  L^^,  and  L^,  are 
much  bigger  than  the  other  elements  in  the  same  row;  most  of  the  non¬ 
diagonal  elements  are  zero.  The  submatrices  are  of  equal  size. 

The  matrix  of  coefficients  L  is  of  the  form  (5.2b) 
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The  diagonal  and  non-diagonal  elements  of  each  part  can  be 


calculated  as  follows, 

i^pq  -  e°pGpq  -  f°pBpq  (5.3) 

l^pp  =  OpGpp  -  fp3pp  +  CIRp  (5.4) 

lI3:pq  =  e°PBPP  +  f°PGpq  (5.5) 

l^^pp  =  epBpp  +  -^pGpp  G^-p  (5.6) 


1 

1 


III 

pq 


hi 

pp 


e°pBpq  r  f°pGpq 

epBpp  +  fpGpp  -  GUp 


lIVpq  =  £°pBpq 


pO  r 

e  p1-  pq 


(5.7) 

(5.8) 


(5.9) 
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(5.10) 


There  are  some  equalities  that  help  in  deciding  whether  the 
condition  (5.1)  is  satisfied  or  not. 

From  the  computation  technique  of  the  bus  admittance  matrix, 
(see  Appendix  2) ,  it  is  evident  that 


(5.11) 


(5.12) 


N 


In  (5.12)  b’pq  is  the  line  charging  to  ground,  and  Cp 
would  be  a  synchronous  condenser  connected  to  bus  "p".  These 


relations  can  be  altered  by  the  presence  of  a  variable  tap 
setting  transformer  between  node  "p",  and  any  other  node. 

From  the  nature  of  the  expressions  (5.3)  to  (5.10),  the 
range  of  values  and  signs  that  CIRp  and  Clip  may  have,  and  the 
fact  that  e°[,  and  f°j_  are  supposed  to  be  close  to  the  exact 
values  eqq,  and  f*q,  it  can  be  concluded  that,  as  far  as  each 

submatrix  is  concerned,  the  diagonal  elements  will  have  an 
absolute  value  close  to  the  absolute  value  of  the  sum  of  the 
other  terms  in  the  same  row. 

In  case  that  a  bus  "p"  is  connected  to  a  slack  bus,  the 
diagonal  term  will  be  predominant,  i.e.:  its  absolute  value 
will  be  bigger  than  the  absolute  value  of  the  sum  of  the  terms 
in  the  same  row.  This  can  be  explained  by  the  fact  that  the 
equations  (5.11)  and  (5.12)  are  always  satisfied,  the  terms 
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that  would  be  called  lpS,  ("s"  standing  for  slack  bus),  are 
missing  because  they  are  zero. 

Hence,  it  can  be  seen  that  each  submatrix  L^,  L^,  L111, 
LIV,  complies  with  the  condition  (5.1),  and  for  the  rows  that 
correspond  to  a  bus  directly  connected  to  a  slack  bus,  the 
condition  (5.1)  is  clearly  satisfied.  Hox^ever,  the  problem 
deals  with  L  as  a  whole,  therefore  we  shall  see  how  L  can 
satisfy  the  condition  (5.1)  as  well. 

Each  one  of  the  submatrices  L^,  L11,  and  L^, 

has  predominant  diagonal  elements.  For  the  application  of 
the  ordinary  Gauss-Seidel  technique  the  convergence  is 
guaranteed  if  the  main  diagonal  elements  of  L  are  dominant. 

For  a  bus  "p"  it  could  be  that  lpp  and  lx  pp,  which  have 
similar  value,  are  more  predominant  that  l^pp  and  l-*--^  ,  whose 

numerical  values  are  close  to  each  other.  Although  condition 
(5.1)  is  not  necessary  for  the  convergence  of  the  Gauss-Seidel 
technique,  applying  a  procedure  that  is  capable  of  transform¬ 
ing  the  linearized  system  in  such  a  way  that  condition  (5.1) 
is  satisfied,  will  bring  a.bout  better  numerical  results. 

5. 1.2.1  The  elimination  of  the  diagonal  of  The 

satisfaction  of  the  condition  (5.1)  could  be  improved  in  several 
ways  if  the  diagonal  terms  of  L11  and  L111  were  to  be  more 
dominant  than  those  of  L1  and  One  of  the  ways  would  be 

by  the  rotation  of  the  submatrices,  puting  L  at  the  top 
left  side  and  L1  at  the  top  right  side,  and  then  LIV  at  the 
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bottom  left  side  and  L"*'^  at  the  bottom  right  side.  The  terms 
of  the  vector  of  increment  of  voltage,  that  come  from  the 
solution  of  the  linearized  systems,  would  then  be  in  a  differ¬ 
ent  order:  Afl*  Af2>***>  AfN>  Aeq,  Ae2,...,  Ae^.  This  method 
has  not  been  tried  because  other  less  complicated  methods,  as 
far  as  computer  programming  is  concerned,  have  proved  to  work 
well  enough. 

The  Gauss-Seidel  method  is  a  succesive  displacement  method, 
which  means  that  the  recently  calculated  elements  of  the 
solution  vector  are  successively  used  as  they  are  obtained. 

Since  the  diagonal  elements  are  predominant  in  all  four 
submatrices,  it  is  our  purpose  to  find  an  algorithm  such  that  a 
simultaneous  solution  for  Aeq  and  Afq  could  be  found  with  it.  In 
this  way  the  error  that  would  be  brought  about  by  multiplying 
the  big  diagonal  elements  by  the  inexact  assumed  values  coming 
from  the  former  run,  would  be  lessened. 

The  algorithm  used  for  the  implementation  of  the  Gauss- 

Seidel  technique,  and  for  approaching  the  linearized  system  to 

the  satisfaction  of  the  condition  (5.1),  consists  of  the 

ttt 

successive  elimination  of  the  diagonal  elements  of  L1  .  The 
elimination  is  nothing  but  a  set  of  equivalent  transformations 
of  the  original  linearized  system. 

For  illustration,  the  operations  will  be  made  on  the  matrix 
pattern  (5.2b).  At  first,  the  first  diagonal  element  of  K, 

will  be  eliminated  by  substracting  from  its  row,  the  first  row  of 
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the  matrix  L,  where  the  diagonal  term  A  lies,  multiplied  by  the  factor 
K/A.  Then  the  increment  Afp  is  calculated  by  multiplying  all  the 
elements  of  the  transformed  row  except  (L  -  KxB/L) ,  which  is  the  trans¬ 
form  of  element  L,  by  the  estimations  Aep,  and  Afp,  on  a  ordered,  ele¬ 
ment  by  element  basis.  Then  the  sum  of  all  these  products  is  subs- 
tracted  from  the  transformed  AQp.  Dividing  the  subtraction  by 
(L  -  KxB/L),  a  new  value  for  Afp  is  obtained.  So  far  no  diagonal  ele¬ 
ment  has  been  multiplied  by  any  Aep  or  Afp,  because  K  was  eliminated. 

Aep  is  calculated  by  multiplying  ^  >  ag ,  a^,  a^ ,  B,  b2?  bg,  b^, 
b5,  by  Aegj  Ae^ ,  Ae^,  Ae^,  Afp,  Af2,  Af^,  Af^,  Af^,  respectively.  The 
products  are  then  added  up,  and  the  addition  substracted  from  APp. 

The  result  is  then  divided  by  A  thus  obtaining  the  new  value  Aep.  It 
can  be  noticed  that  in  the  calculation  of  Aep,  one  big  diagonal  ele¬ 
ment,  B,  has  been  multiplied,  but  it  has  been  multiplied  by  Afp,  which 
was  calculated  just  at  the  step  before,  therefore  being  closer  to  the 
solution  than  the  previous  value. 

The  same  is  now  done  with  the  second  diagonal  element  of 
M  is  eliminated  by  substracting  the  second  row  of  the  matrix  multi¬ 
plied  by  M/C;  Afp  is  calculated,  and  then  a  new  Ae2  is  computed  using 
already  the  values  of  Af^,  and  Aep  and  Afp,  obtained  in  the  previous 
step.  The  same  process  is  repeated  until  a  whole  new  set  Aep,  Afp 
has  been  calculated,  and  so  completed  a  run  of  the  Gauss-Seidel  tech¬ 
nique  . 

Since  the  constant  elements  AQp,  AQ2,...,  Aq^,  are  modified 
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during  the  elimination  of  the  elements  of  the  diagonal  of  L111, 
and  they  are  long  to  compute,  an  extra  vector  to  keep  those  values 
must  be  supplied. 

5. 1.2. 2  Characteristics  of  the  method.  The  Gauss-Seidel 
technique,  using  the  algorithm  described  above,  has  some  ad¬ 
vantages  . 

-  The  algorithm  is  very  easy  to  program,  because  it  is  only 
a  sequence  of  operations,  which  can  be  solved  in  a  single  loop. 

-  In  case  that  the  diagonal  terms  l^pp  and  l^^p  are 

more  dominant  that  l^p  and  l^pp,  the  elimination  of  the 
diagonal  of  yields  an  equivalent  system.  In  this  system  the 

new  l^pp  elements  are  more  predominant  than  before,  because  the 
new  l^pp  is  a  linear  combination  of  the  ancient  1-^p  and  l-^pp, 
which  is  more  diagonal  dominant.  It  can  be  said  that,  since  the 
rows  of  the  lower  half  of  L  are  transformed  into  a  row  with  only 
one  dominant  element,  the  condition  (5.1)  will  be  satisfied  more 
largely . 

-  The  more  slack  buses  are  used,  the  faster  the  convergence 
of  the  Gauss-Seidel  technique.  As  explained  before,  the  buses 
connected  to  a  slack  bus  have  the  diagonal  term  predominant, 
because  in  its  row,  the  term  corresponding  to  the  slack  bus  is 
missing.  This  is  an  advantage  over  the  triangulation  with 
pivoting  technique,  where  the  number  of  slack  buses  does  not 
affect  the  speed  of  solution. 

-  Another  advantage  over  the  triangulation  with  pivoting 
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technique  is  that,  as  stated  in  section  4  of  chapter  2  in  refer¬ 
ence  (11) ,  errors  due  to  roundoff  and  blunders  may  be  damped  out 
as  the  Gauss-Seidel  technique  proceeds. 

5.2  THE  APPLICATION  OF  AN  ITERATIVE  TECHNIQUE  FOR  SOLVING  THE 

LINEARIZED  SYSTEM. 

Experimental  results  with  the  IEEE  14,  30,  and  57  bus  test 
systems  show  that  the  Gauss-Seidel  technique  is  convergent  but 
its  speed  is  very  low  .  Acceleration  factors  have  proved  to 
increase  the  speed  of  the  convergence  of  the  Gauss-Seidel  tech¬ 
nique  for  solving  the  linearized  system.  Even  though  the  speed 
of  convergence  is  low,  the  proposed  method  with  the  Gauss- 
Seidel  technique  gives  faster  results  than  the  Newton-Raphson 
method  using  triangulation  with  pivoting,  in  some  cases. 

5.2.1  DIFFERENCE  BETWEEN  ITERATION  AND  RUN.  In  order  to  avoid 
confusion  a  different  name  has  been  used  for  the  iterations  in 
the  proposed  method  or  in  the  Newton-Raphson  method,  (which  have 
been  called  "iterations"  all  along) ,  and  for  the  iterations  in 
the  Gauss-Seidel  technique  for  solving  the  linearized  system. 
Each  of  the  iterations,  namely,  the  calculation  of  a  new  set  of 
values  Aei>  Afi>'  has  been  called  a  "run".  For  instance,  we 
could  use  10  "runs"  of  the  Gauss-Seidel  technique  at  each 
iteration,  which  would  mean  that  the  process  of  obtaining  a  new 
set  of  values  Aej_,  Afp,  would  be  applied  ten  times  between  two 
changes  of  the  diagonal  elements  of  the  coefficient  matrix. 
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The  number  of  "runs"  can  be  constant,  or  variable,  or 
regulated  by  a  criteria,  which  measures  the  rate  of  convergence 
of  the  results  in  the  runs,  and  stops  the  calculations  in  a  new 
run  when  the  rates  of  convergence  reach  a  specified  value. 

When  initiating  an  "iteration",  a  check  of  the  convergence 
in  voltage  and  in  power  is  made,  and  the  values  of  the  voltage 
and  the  currents  in  the  diagonals  of  the  linearized  system  matrix 
L  are  updated.  Once  this  is  done,  the  program  goes  to  the  routine 
for  solving  the  new  linear  system.  If  this  routine  is  a  Gauss-Seidel 
technique  such  as  described  above,  a  certain  number  of  "runs"  will 
be  used  to  get  a  new,  better  approximation  of  the  solution. 

When  the  Gauss-Seidel  technique  is  used  for  solving  the  system 
of  linear  equations,  a  starting  set  of  values  must  be  supplied.  In 
the  proposed  method  the  set  of  starting  values  is  only  supplied  at 
the  first  iteration.  These  values  are  Ae-j_  =  0.0  ,  Af-j_  =  0.0  , 
because  it  is  not  known  in  which  direction,  with  respect  to  the 
initial  values  e0-^,  f°i,  the  solution  is  going  to  be.  During  the 
iterations  other  than  the  first,  the  result  of  the  former  iteration 
is  the  best  starting  value  that  can  be  used. 

If  the  Gauss-Seidel  technique  along  with  the  Newton-Raphson 
method  had  been  used,  the  starting  values  would  have  to  be  zero  at 
the  beginning  of  each  iteration,  because  the  increments  at  the  be¬ 
ginning  of  a  new  iteration  have  nothing  to  do  with  the  increments 
of  voltage  in  the  former  iteration. 
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obtained.  Applying  this  process  a  few  times,  the  precision  obtained 
in  the  first  iteration  using  either  triangulation  with  pivoting  or 
a  high  number  of  runs  per  iteration  with  the  Gauss-Seidel  technique, 
is  reached,  and  the  time  elapsed  proves  to  be  less. 

Another  advantage  is  the  possible  application  of  acceleration 
factors,  which  can  be  implemented  at  any  stage  of  the  Gauss-Seidel 
technique.  In  chapter  7  is  shown  that  the  acceleration  factors 
speed  up  the  solution. 
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5.2.2  TENDENCY  OF  THE  RESULTS  OBTAINED  IN  THE  RUNS.  As  it  was 
pointed  out  before,  it  is  a  waste  of  time  to  solve  the  linearized 
system  with  high  precision  at  the  first  iteration  because  the 
values  CIRp  and  Clip  are  far  from  the  final  values.  Time  would  be 
spent  for  solving  with  accuracy  an  approximate  system.  Even  being 
solved  exactly,  it  would  only  give  the  tendency  towards  the  solution. 

It  has  been  found  that  the  same  tendency  towards  the  solution 
can  be  obtained  in  less  time  using  the  Gauss-Seidel  technique  for 
solving  the  linearized  system,  but  using  only  a  few  runs  instead  of 
a  complete  solution,  which  would  take  a  lot  of  runs.  The  complete 
solution  by  the  Gauss-Seidel  technique  takes  a  much  longer  time 
than  the  triangulation  with  pivoting  technique,  due  to  the  low  speed 
of  convergence. 

Advantage  can  be  taken  of  two  facts:  first,  whereas  the  trian¬ 
gulation  technique  has  to  be  carried  out  completely  no  matter  how 
much  accuracy  is  needed,  with  the  Gauss-Seidel  technique,  the  cal¬ 
culation  of  "runs"  can  be  stopped  at  any  time,  so,  if  no  high 
precision  is  needed,  a  few  runs  will  be  enough.  Of  course  the  vol¬ 
tages  are  not  as  close  to  the  solution  after  a  small  number  of  runs 
as  they  would  have  been  if  more  runs  per  iteration  or  a  triangula¬ 
tion  routine  had  been  used.  Since  the  tendency  of  the  results  of  the 
runs  is  always  towards  the  solution,  as  the  system  is  convergent,  if 
the  currents  CIRp  and  Clip  are  recomputed,  better  values  for  them 
would  be  obtained.  Substituting  them  in  the  diagonal  elements  of  the 
linearized  system,  a  system  that  will  yield  an  improved  solution  is 
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VI. 


RESULTS  OF  THE  PROPOSED  METHOD  USING  THE  GAUSS- 


SEIDEL  ITERATIVE  TECHNIQUE  FOR  SOLVING  THE 

LINEARIZED  SYSTEM 

6.1  ANALYSIS  OF  THE  RESULTS  OF  AN  EXACT,  AND  AN  APPROXIMATE 

SOLUTION  OF  THE  LINEARIZED  SYSTEM  IN  THE  PROPOSED  METHOD 

In  section  5.2.2  it  was  reported  that  an  approximate  iterative 
solution  of  the  linearized  system  would  yield  time  advantages.  The 
approximate  solution  is  obtained  by  using  only  a  few  runs  of  the 
Gauss-Seidel  technique  instead  of  as  many  as  necessary  for  the  exact 
solution  of  the  linearized  system  at  each  iteration. 

For  the  solution  of  the  linearized  system  at  each  iteration  a 
different  number  of  runs  could  be  used,  but  only  the  results  with  a 
constant  number  of  runs  per  iteration  will  be  shown.  One  of  the 
improvements  of  the  proposed  method  could  be  the  increasing  of  the 
number  of  runs  per  iteration  at  each  iteration,  so  as  to  obtain  more 
accurate  results  as  the  solution  is  approached.  However,  even  a 
constant  number  of  runs  per  iteration  improves  the  solution  time. 

All  the  results  use  the  same  degree  of  precision  so  that  a 
fair  comparison  can  be  drawn.  It  must  be  pointed  out  that 
when  using  the  IBM  library  subroutines  CS019A  and  CS019B  for 
measuring  the  time  spent,  account  has  to  be  taken  of  several  facts: 
first,  the  amount  of  output.  All  printing  operations  are  time 
consuming,  and  the  imput-output  modes  use  a  variable  time  depending 
on  the  circumstances.  To  avoid  that  each  case  has  been  solved 
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twice,  one  time  with  printing  of  the  intermediate  results,  and 
another  without  any  printing  until  the  final  results  are  reached. 

Even  so,  sometimes  the  time  employed  for  the  same  problem  is 
different  due  to  the  fact  that  the  computer  works  on  a  shared  time 
basis.  To  overcome  that  there  has  been  run  along  with  all  the 
cases  a  reference  case,  which  consists  of  the  Newton-Raphson 
solution. 

The  proposed  method  for  solving  the  load-flow  problem  is 
convergent  for  the  IEEE  14,  30  and  57  bus  test  systems,  and, 
besides,  as  specified  in  section  5.12,  the  iterative  technique  for 
the  solution  of  the  linearized  system  is  convergent  too.  Several 
cases  using  different  runs  per  iteration  are  shown  for  each  test 
system. 

Since  the  results  of  the  load-flow  problem  for  each  system  is 
known,  the  convergence  at  each  step  of  the  solution  can  be  checked 
by  finding  the  differences  between  the  final  voltages  and  the  vol¬ 
tages  at  each  stage  of  the  solution. 

6.1.1  RESULTS  FOR  THE  14  BUS  TEST  SYSTEM. 

In  table  6.1  is  shown  that  there  exists  a  certain  number  of  runs 
per  iteration  which  gives  the  solution  in  minimum  time  with  respect 
to  other  numbers  of  runs  per  iteration.  This  is  so  because  as  the 
number  of  runs  per  iteration  decreases,  the  time  per  iteration 
decreases,  but  the  number  of  iterations  for  convergence  increases. 

The  minimum  is  the  best  compromise. 
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TABLE  6.1  TIMES  AND  NUMBER  OF  ITERATIONS  FOR  SEVERAL  CASES  OF 
APPROXIMATE  SOLUTION  OF  THE  LINEARIZED  SYSTEM  IN  THE 
PROPOSED  METHOD  APPLIED  TO  THE  14  BUS  SYSTEM. 


runs  per 
iteration 

iterations 
for  convergences 

time 

solution 

150 

6 

10.7 

40 

6 

4.54 

20 

6 

2.237 

13 

5 

1.488 

10 

7 

1.641 

8 

10 

1.918 

6 

12 

1.801 

5 

14 

1.811 

4 

17 

1.97 

3 

22 

2.144 

2 

30 

2.132 

1 

50 

2.3 

Newton-Raphs on 

4 

0.624 

for 

(seconds) 


N.B.  The  program  used  deals  with  the  whole  Y^us  matrix,  i.e.: 
there  are  no  logical  instructions  in  the  program  to  deal  with 
only  the  non-zero  elements  of  Y^us. 


Some  of  the  maximum  and  minimum  active  and  reactive  voltage 
error  with  respect  to  the  final  values  are  shown  in  table  6.2. 

It  can  be  observed  that  the  error  at  the  end  of  each  iteration 
gets  bigger  as  the  number  of  runs  per  iteration  used  decreases. 
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TABLE  6.2  MAXIMUM  AND  MINIMUM  VOLTAGE  ERRORS  AT  EACH  ITERATION  OF  THE  PROPOSED  METHOD  WITH  SEVERAL 
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The  results  in  table  6.  1  reveal  that  even  with  only  one  run  per 
iteration  the  proposed  method  is  convergent.  The  convergence  speed 
decreases  as  predicted  according  to  the  theoretical  results  in 
sections  3.1.1  and  3.3. 

The  time  spent  using  10  and  4  runs  per  iteration  .  is  shown  in 
table  6.3  .  The  program  used  deals  only  with  the  nonzero  elements 
of  the  matrix  Ybus,  therefore  the  program  has  logical  instructions 
to  work  only  with  those  elements.  The  results  could  be  compared 
with  those  in  table  2.1,  which  correspond  to  the  Newton-Raphson 
method . 


TABLE  6.3  TIME  SPENT  BY  THE  PROPOSED  METHOD  USING  10  AND  4  RUNS 
PER  ITERATION  WHEN  APPLIED  TO  THE  14  BUS  TEST  SYSTEM 


10  runs  per  iteration,  times  in  seconds 


iteration 

number 


time  per  new  currents  new  terms  of 
iteration  and  powers  linearized  sys . 


sys .  solution 
and  new  volts. 


1  0.3927  0.0175  0.1613 


0.2099 


2-7  0.2312  0.0169  0.0028 


0.2115 


4  runs  per  iteration,  times  in  seconds 


iteration 

number 


time  per  new  currents  new  terms  of 
iteration  and  powers  linearized  sys. 


sys.  solution 
and  new  volts. 


0.2656 


0.0177 


0.1613 


0.0866 


2-17 


0.1074 


0.0171 


0.0041 


0.0863 


N.B.  The  times  spent  for  the  calculation  of  the  new  currents  and 
powers,  and  for  the  solution  of  the  linearized  system  and 
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calculation  of  the  new  voltages  ought  to  be  the  same  at  each 
iteration  for  each  case;  however,  small  variations  are  observed. 
Moreover,  small  variations  are  observed  as  well  in  the  time  spent 
on  the  different  items  of  the  iterations  beyond  the  first  one. 

Comparing  the  times  taken  at  each  step  of  the  iteration,  as 
shown  in  table  6.3,  with  the  times  taken  by  the  Newton-Raphson  me¬ 
thod,  the  following  can  be  concluded: 

-  The  time  for  computing  the  new  currents  and  powers  is  more 
or  less  the  same  using  the  Newton-Raphson  method  and  the  proposed 
method,  and  in  the  last,  it  is  independent  of  the  number  of  runs  per 
iteration,  as  expected. 

-  In  the  proposed  method,  the  time  for  the  computing  of  the 
new  elements  of  the  linearized  system  is  long  in  the  first  itera¬ 
tion  and  short  in  the  others.  The  fact  of  calculating  only  the 
diagonal  elements  after  the  first  iteracion  accounts  for  it.  The 
time  spent  on  this  step  in  the  first  iteration  of  the  proposed 
method  is  more  or  less  the  same  as  the  corresponding  time  for  the 
Newton-Raphson  method.  Here  lies  one  of  the  advantages  of  the 
proposed  method:  this  method  does  not  spend  time  in  recalculating 
all  the  coefficients  of  the  linearized  system  at  each  iteration. 

-  The  time  for  solving  the  linear  system  is  variable  depend¬ 
ing  on  the  number  of  runs  per  iteration:  0.2099  s.  for  10  runs  per 
iteration,  and  0.0866  s.  for  4  runs  per  iteration.  The  triangula¬ 
tion  with  pivoting  used  with  the  Newton-Raphson  method  takes 
0.1662  s.  .  Therefore,  depending  on  the  number  of  runs  per 
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iteration,  time  will  be  saved  or  lost  with  respect  to  the  Newton- 
Raphson  method. 

-  A  comment  on  the  precision  can  be  found  with  the  57  bus 

case. 


6.1.2  RESULTS  FOR  THE  30  BUS  TEST  SYSTEM 

TABLE  6.4  TIMES  AND  NUMBER  OF  ITERATIONS  FOR  SEVERAL  CASES  OF 

APPROXIMATE  SOLUTIONS  OF  THE  LINEARIZED  SYSTEM  IN  THE 
PROPOSED  METHOD  APPLIED  TO  THE  30  BUS  SYSTEM. 


runs  per 
iteration 


iterations  time  for 

for  converg.  solution  (seconds) 


150 

6 

— 

40 

7 

- 

20 

9 

- 

13 

11 

14.09 

10 

15 

15.03 

8 

18 

15.08 

6 

23 

14.70 

5 

26 

14.49 

4 

31 

14.16 

3 

38 

14.08 

2 

49 

13 . 64 

1 

73 

13.83 

N-R 


4 


4.45 


The  same  N.B.  as  in  Table  6.1  applies  here. 

Table  6.5  shows  the  voltage  errors  for  this  test  system. 
The  times  obtained  when  using  the  proposed  method  with  5 
runs  per  iteration  are  shown  in  table  6.6  . 
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TABLE  6.6  TIME  SPENT  BY  THE  PROPOSED  METHOD  USING  5  RUNS  PER 


ITERATION  WHEN  APPLIED  TO  THE  30  BUS  TEST  SYSTEM. 


Times  in  seconds 


iteration 

number 


time  per 
iteration 


new  currents  new  terms  of 
and  powers  linearized  sys . 


sys .  solution 
and  new  volts. 


1  2.0566  0.0527  1.5603 

2-26  0.5028  0.0518  0.0057 


0.4436 

0.4453 


The  same  N.B.  as  in  table  6.3  applies  here. 

Comparing  the  results  in  Table  6.6  with  the  results  of  the 
Newton-Raphson  method  for  the  same  system,  the  same  conclusions  as 
for  the  14  bus  test  system  would  be  reached. 


6.1.3  RESULTS  FOR  THE  57  BUS  TEST  SYSTEM 

TABLE  6.7  TIMES  AND  NUMBER  OF  ITERATIONS  FOR  SEVERAL  CASES  OF 

APPROXIMATE  SOLUTIONS  OF  THE  LINEARIZED  SYSTEM  IN  THE 
PROPOSED  METHOD  APPLIED  TO  THE  57  BUS  SYSTEM. 


runs  per 
iteration 

iterations 
for  converg. 

time  for 
solution 

time  for 
first  ite. 

time  for 
iteration 

24 

14 

122.8 

19.74 

8.00 

18 

18 

125.2 

17.9 

6.16 

12 

25 

121.0 

15.85 

4.16 

10 

28 

106.0 

14.9 

3.35 

8 

34 

111.0 

14.54 

2.9 

6 

42 

103.0 

13.83 

2.165 

5 

48 

102.9 

13.48 

1.9 

4 

56 

96.3 

13.2 

1.5 

. 
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runs  per 
iteration 


iterations  time  for 
for  converg.  solution 


time  for  time  for 

first  ite.  iteration  (seconds) 


2 

N-R 


12.2 

4  (3  compl.)  60.0  19.88 


0.835 

19.88 


The  program  used  deals  only  with  the  non-zero  elements  of  Ybus 
by  means  of  logical  instructions. 

It  can  be  observed  in  table  6.7,  that  the  times  for  the  first 
iteration,  are  much  higher  than  those  for  the  other  iterations.  The 
reason  is  that  during  the  first  iteration  all  the  elements  of  the 
linearized  system  must  be  calculated. 

The  voltage  errors  for  the  57  bus  test  system  are  found  in 
table  6.8  . 


TABLE  6.8  MAXIMUM  AND  MINIMUM  VOLTAGE  ERRORS  AT  EACH  ITERATION  OF 
THE  PROPOSED  METHOD  WITH  SEVERAL  CASES  OF  APPROXIMATE 
SOLUTIONS  OF  THE  LINEARIZED  SYSTEM.  APPLICATION  TO  THE 
57  BUS  TEST  SYSTEM. 

The  same  N.B.  as  in  table  6.2  applies  here. 


runs 

per  ite 

.  1 

ite 

.  2 

ite 

.  3 

ite . 

4 

ite .  5 

iteration 

-1301 

-803 

-63 

-29 

0 

0 

N-R 

-4782  - 

2836 

-212 

-142 

0 

0 

935 

9348 

661 

7380 

460 

5664 

333 

4784 

10 

-6770 

458 

-4511 

-336 

-3280 

460 

-2647 

378 

1196 

9876 

863 

8214 

618 

6181 

464 

5488 

355  4803 

8 

-7357 

670 

-5243 

-173 

-3855 

584 

-3026 

488 

-2532  353 

' 


‘ 
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runs  per 
iteration 

ite.  1 

ite. 

4 

ite. 

8 

ite 

.  16 

149  7 

10528 

662 

6527 

6 

7993 

969 

-3787 

646 

1822 

11393 

989 

8104 

518 

5512 

4 

■8680 

651 

-5165 

888 

-2826 

477 

2468 

12551 

1599 

9777 

1037 

8098 

541 

5532 

2 

■9412 

208 

-7315 

1373 

-5120 

877 

-2727 

487 

TABLE  6.9 

TIME 

SPENT 

BY  THE  PROPOSED 

METHOD 

USING 

12  RUNS 

PER 

ITERATION  WHEN  APPLIED  TO  THE  57  BUS  TEST  SYSTEM. 


Times  in 

seconds 

iteration 

number 

time  per 
iteration 

new  currents 
and  powers 

new  terms  of 
linearized  sys. 

sys.  solution 
and  new  volts 

1 

15.2760 

0.1448 

11.2762 

3.8550 

2-25 

4.0192 

0.1434 

0.0110 

3.8648 

The  same  N.B.  as  in  table  6.3  applies  here. 

6.2  COMPARISON  BETWEEN  THE  RESULTS  IN  THE  THREE  SYSTEMS  AND 
CONCLUSIONS 

The  following  facts  can  be  observed: 

-  The  values  of  the  errors  at  each  iteration  for  the  57  bus 
case  are  much  higher  than  for  the  14  and  30  bus  cases.  This  is  due 
to  several  reasons.  At  first,  the  initial  guess  E^  =  1.0,  Fjc  =  -0.3 
is  not  as  good  in  this  case  as  it  is  for  the  14  and  30  bus  case. 

It  is  very  likely  that  F^  =  -0.22  would  have  worked  better,  for 
this  value  is  closer  to  the  average  final  value  for  the  reactive 
voltage.  Moreover  the  index  number  of  lines/number  of  nodes  is 
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1.403,  (1.366  for  the  30  bus  system),  which  affects  the  error  at  the 
end  of  the  first  iteration.  Indeed,  during  the  first  run  of  the  first 
iteration  of  the  proposed  method,  a  starting  value  of  zero  for  all  the 
increments  is  used.  Therefore,  the  more  lines  conected  to 
each  bus,  the  bigger  the  relative  error  is  at  the  end  of  the  first 
run. 

-  The  rate  of  convergence  of  the  proposed  method  is  decreasing, 
and  at  the  end  is  much  slower  than  that  of  the  Newton-Raphson  method. 
However,  since  the  time  spent  for  iteration  is  smaller,  there  is  a 
period  of  time  in  which  the  proposed  method  gives  better  results 
than  the  Newton-Raphson  method.  See  section  3.1.3. 

-  A  consequence  of  the  decreasing  convergence  is  that  much 

time  is  spent  on  the  last  iterations,  but  there  is  not  much  improvement 
of  the  solution.  The  program  cuts  off  the  problem  when  the  changes 
in  voltage  are  less  than  the  precision  used,  i.e.i  0.0001.  This  is 
done  in  order  not  to  spend  time  on  an  improvement,  which  is  very 
slowT  and  is  likely  to  be  overcome  in  part  by  round-off. 

-  Another  consequence  of  the  decreasing  convergence  is  that 
although  the  increments  of  voltage  are  less  than  the  tolerance 
0,0001  when  the  solution  is  cut  off,  the  precision  in  power  might 
not  be  below  the  tolerance.  The  accumulated  numerical  error  in  the 

diagonal  terms  of  the  linearized  system,  due  to  a  high  number  of 
iterations,  could  account  for  it. 

A  display  of  the  maximum  differences  in  power  and  the  maximum 
changes  of  voltage  at  the  last  iterations,  can  be  found  in  table  6.10. 


TABLE  6.10  MAXIMUM  DIFFERENCES  IN  POWER  AND  BIGGEST  VOLTAGE  CHANGE  AT  THE  FINAL  ITERATIONS 


79 


CO 

d 

P 

t— 

ir\ 


w 

• 

P 

H 

CU 

P 

O 

O 

cu 

w 

d 

co 

•H 

<3 

u 

CD 

O 

• 

d 

w 

0) 

H 

p 

M 

d) 

4-1 

m 

00 

00 

CM 

ON 

Pi 

CO 

4-1 

CnI 

O 

1 — 1 

rH 

O 

cn 

p 

•H 

o 

o 

o 

o 

p 

P 

• 

o 

o 

• 

o 

o 

P 

d) 

o 

o 

0) 

o 

o 

B 

p 

• 

• 

P 

• 

• 

MO 

d 

■H 

o 

o 

•H 

o 

o 

6 

w 

•H 

3 

H 

d 

<fr 

O 

CO 

rH 

m 

CM 

z 

cu 

CM 

i—r 

CM 

rH 

i — i 

M0 

<2 

P 

o 

O 

o 

O 

4-1 

• 

o 

O 

• 

o 

O 

s 

d 

o 

O 

d) 

o 

o 

w 

P 

4-J 

• 

• 

P 

• 

• 

H 

o 

•H 

o 

o 

•H 

o 

o 

CO 

4-1 

>-l 

CO 

H 

CO 

CO 

w 

CO 

CO 

CO 

CO 

O 

co 

r-. 

H 

d 

CM 

rH 

CO 

i — 1 

CM 

C T\ 

d 

O 

o 

o 

O 

CO 

• 

o 

o 

• 

o 

o 

p 

0) 

d) 

O 

o 

CU 

o 

o 

w 

00 

-U 

• 

• 

p 

• 

• 

d 

•H 

o 

o 

•H 

o 

o 

4-J 

in 

rH 

o 

w 

> 

w 

£h 

d 

•H 

CM 

MO 

CO 

P 

CM 

i — i 

CTv 

CO 

m 

O 

CD 

o 

o 

o 

i — i 

00 

• 

o 

o 

• 

o 

o 

w 

d 

d> 

o 

o 

<U 

o 

o 

CO 

d 

4J 

• 

• 

P 

• 

• 

<3 

P 

•H 

o 

o 

•H 

o 

o 

O 

o 

• 

• 

B 

• 

W 

£ 

d 

a) 

H 

w 

s 

4-J 

H 

Eh 

•rH 

•H 

CO 

X 

i — 1 

1 — 1 

M0 

' — . 

r^ 

I — 1 

CO 

>H 

d 

CM 

CM 

M0 

CO 

CXI 

in 

<f 

1 — ^ 
g 

CO 

B 

• 

O 

o 

d 

o 

CN 

CU 

• 

O 

o 

d 

• 

o 

O 

p 

EH 

<U 

4-1 

d) 

O 

o 

p 

d) 

o 

O 

CO 

P 

*rl 

4J 

• 

• 

P 

• 

• 

CM 

W 

4-) 

•H 

o 

o 

M0 

•H 

o 

o 

rH 

H 

CO 

P 

d 

W 

CO 

o 

d 

B 

P 

p 

4H 

P 

a) 

H 

(O 

p 

CO 

CM 

CO 

P 

<* 

P 

i — 1 

O 

vO 

I — 1 

Jo 

o 

o 

O 

rH 

d 

CM 

CM 

CO\ 

CO 

O'! 

CO 

d 

o 

o 

o 

CO 

4-J 

B 

• 

o 

o 

p 

• 

o 

o 

CO 

0) 

d) 

o 

o 

CO 

CU 

o 

o 

4-1 

P 

. 

• 

d) 

p 

• 

• 

CO 

•rl 

o 

o 

p 

•H 

o 

o 

Jg 

P 

H 

CO 

CO 

d 

4-J 

p 

CO 

X 

X 

£ 

a) 

2 

y 

4-J 

Eh 

CO 

1 — 1 

fcH 

CO 

*  '1 

> 

-80 


From  the  results  in  table  6.10  the  following  can  be  concluded: 

a)  The  convergence  is  indeed  slower  at  the  end  of  the  process. 

b)  There  is  a  trend  towards  a  bigger  mismatch  of  power  as  the  number 
of  runs  per  iteration  decreases.  This  phenomenon  is  general. 

c)  Since  TMX  is  measured  after  the  corrections  of  voltage  and  SMX 
is  measured  before,  the  mismatch  is  not  as  big  as  it  appears 
because  the  value  of  SMX  at  the  iteration  "n+l"  would  in  fact 
correspond  to  the  value  of  TMX  at  the  iteration  "n"  ,  and  SMX 
decreases . 

d)  As  a  consequence  of  a) ,  b) ,  and  c) ,  the  convenience  of  finishing 
off  with  one  iteration  using  the  Newton-Raphson  method,  can  be 
seen.  This  iteration  would  eliminate  the  power  mismatch,  and 
would  converge  more  rapidly. 


VII. 


IMPROVEMENTS  IN  THE  PROPOSED  METHOD 


7.1  THE  USE  OF  ACCELERATION  FACTORS  IN  THE  SOLUTION  OF  THE 
LINEARIZED  SYSTEM  BY  THE  GAUSS-SEIDEL  TECHNIQUE 
Acceleration  is  only  implemented  in  the  solution  of  the 

linearized  system  at  each  iteration,  so  as  to  obtain  more  exact 
results  in  less  iterations.  By  doing  so,  the  extra  time  that 
would  have  been  taken  by  the  non-accelerated  iterations  to  reach 
the  same  precision,  is  saved. 

Other  acceleration  procedures  could  be  used  between  the  values 
at  the  end  of  each  iteration,  but  these  are  the  ordinary  techniques 
that  are  used  with  every  direct  load-flow  method,  and  they  will  not 
be  used  xxrith  either  the  proposed  method  nor  x^ith  the  Nexjton-Raphson 
method,  which  is  taken  as  the  reference  case. 

Some  values  of  acceleration  factors  have  been  tried.  The 
values  used  improve  the  speed  more  than  other  values  tested,  but  it 
is  possible  they  may  be  improved  by  another  set  of  values.  A 
search  for  the  best  values  has  not  been  intended  because  the  speed 
of  the  process  can  still  be  improved  in  many  other  ways,  and  be¬ 
sides,  the  optimum  values  would  be  different  for  each  system.  The 
goal  was  to  show  that  acceleration  factors  improve  the  solution, 
sometimes  to  such  an  extent  that  xve  reach  the  solution  of  the  load 
floxtf  problem  in  less  time  than  the  Newton-Raphson  method. 

7.1.1  ACCELERATION  FORMULA  EMPLOYED.  At  the  end  of  each  run  the 
iterative  technique  comes  up  xvith  a  set  of  increments  of  voltage 
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with  respect  to  the  first  guess.  They  will  be  called  DFT^. .  The  values 
of  DFTk  are  stored  for  three  consecutive  runs,  which  will  be  called: 
r-2,  r-1,  and  r.  Then  a  check  is  made  to  see  if,  for  every  component, 
DFTkr~l  lies  between  DFTkr~2  and  DFTkr .  If  this  is  so,  the  value  of 
the  component  is  modified  according  to  the  formula  (7.1). 

DFTrk  new  =  DFTrk  +  (accl.  f actor)x (DFTkr  -  DFTkr_1)  (7.1) 

A  different  acceleration  factor  is  used  for  the  components  that 
correspond  to  the  active  voltages,  and  those  corresponding  to  the 
reactive  voltages.  Details  about  acceleration  factors  can  be  found 
in  section  2.5  of  reference  (11). 

The  acceleration  formula  (7.1)  is  one  of  those  possible.  Due  to 
the  checking  that  is  made  to  see  whether  DFTkr“^  lies  between  the 
values  at  the  iteration  before  and  after,  the  components  of  DFT 
whose  values  oscillate,  are  not  modified. 

7.2  ANALYSIS  OF  THE  EFFECT  OF  ACCELERATION  FACTORS  IN  THE  SOLUTION 

OF  THE  LINEARIZED  SYSTEM 

The  effect  of  acceleration  factors  is  shown  by  the  performance 
of  the  proposed  method  using  several  values  for  them.  FRL  will  be 
the  acceleration  factor  for  the  real  part  of  the  voltage  increments, 
and  FIM  that  for  the  imaginary  part  of  the  voltage  increments. 

When  the  value  of  the  acceleration  factors  is  FRL  =  0.0,  and  FIM  = 

0.0,  it  is  as  if  no  acceleration  factors  were  applied  to  the  ite¬ 
rative  technique  of  solution  of  the  linearized  system. 
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In  table  7.1  several  pairs  of  acceleration  factors  are  tried 
on  the  3  runs/iteration  case  of  the  solution  of  the  57  bus  test 
system  by  the  proposed  method.  It  can  be  seen  that  for  the  same 
iteration,  the  pair  FRL  3.5  FIM  5.0  gives  the  lowest  error  limits. 

A  graph  of  the  acceleration  factor  effects  is  shown  in  figure 

7.1. 

For  the  30  bus  test  system,  the  sensitivity  with  respect  to 
the  acceleration  factors,  is  shown  by  the  number  of  iterations  and 
time  to  convergence.  In  table  7.2,  the  6  runs  par  iteration  case 
is  shown.  It  can  be  noticed  that  three  of  the  acceleration  factors 
tried,  yield  the  solution  in  less  time  than  the  Newton-Raphson 
method . 
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TABLE  7.2  EFFECT  OF  SEVERAL  ACCELERATIONS  FACTORS  ON  THE  NUMBER  OF 
ITERATIONS  AND  TIME  TO  CONVERGENCE.  APPLICATION  TO  THE 
6  RUNS /ITERATION  CASE  OF  THE  SOLUTION  OF  THE  30  BUS  TEST 
SYSTEM. 


case 

N-R 

FRL  0.0 

4.0 

3.0 

3.0 

3.0 

FIM  0.0 

8.0 

6.0 

6.5 

5.5 

number  of  ite. 

4 

23 

20 

11 

11 

10 

for  convergen. 

time  for  conv. 

8.64 

16.51 

15.51 

8.00 

8.00 

7.71 

(in  seconds) 


In  table  7.3  there  are  the  results  of  several  cases  of  solution 
of  the  14  bus  test  system  by  the  proposed  method  using  acceleration 
factors  for  the  solution  of  the  linearized  system. 

TABLE  7.3  EFFECT  OF  ACCELERATION  FACTORS  ON  THE  TIMES  OF  SOLUTION 
AND  NUMBER  OF  ITERATIONS  FOR  CONVERGENCE.  APPLICATION 
TO  THE  4  RUNS /ITERATION  CASE  OF  THE  SOLUTION  OF  THE  14 
BUS  TEST  SYSTEM. 


Solution  by  the  Newton-Raphson  method  using  triangulation  with 
pivoting  for  the  linearized  system.  Total  time  is  0.9955  s. 


- - 0 - 

:eration 

time  per 

new  currents 

new  terms 

of 

sys.  solution 

number 

iteration 

and  powers 

linearized 

sys . 

and  new  volts. 

1-3 

0.3262 

0.0167 

0.1551 

0.1544 

4 

0.0167 

0.0167 

o 

• 

o 

0.0 
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Proposed  method,  4  runs /iteration.  FRL  =  3.0  ,  FIM  =  6.0  . 


Total  time 

1.4034  s. 

iteration 

number 

time  per 
iteration 

new  currents 
and  powers 

new  terms  of 
linearized  sys . 

system 

solution 

and  new  volts. 

1 

0.2569 

0.0167 

0.1561 

0.0840 

2-12 

0.1041 

0.0163 

0.0031 

0.0847 

Proposed  method,  4  runs 

/iteration.  FRL 

=  3.0  ,  FIM  = 

6.0  . 

Total  time 

0 .9868  s . 

iteration 

number 

time  per 
iteration 

new  currents 
and  powers 

new  terms  of 
linearized  sys 

sys.  solution 
.  and  new  volts 

1 

0.2568 

0.0167 

0.1561 

0.0840 

2-8 

0.1041 

0.0163. 

0.0031 

0.0847 

Proposed  method,  4  runs 

/iteration.  FRL 

=  4.0  ,  FIM  =  8 

.0  . 

Total  time 

1.4029  s. 

iteration 

number 

time  per 
iteration 

new  currents 
and  powers 

new  terms  of 
linearized  sys. 

sys.  solution 
and  new  volts. 

1 

0.2559 

0.0167 

0.1551 

0.0840 

2-12 

0.1041 

0.0163 

0.0031 

0.0847 

The  same  N.B.  as  in  Table  6.3  applies  here. 

7.3  COMMENTS  ON  THE  USE  OF  ACCELERATION  FACTORS 

As  can  be  supposed  from  the  acceleration  formula,  the  subroutine 
which  performs  the  acceleration  works  as  soon  as  three  succesive  sets 


of  DFT  are  given  to  the  subroutine.  One  set  is  given  to  the  sub- 
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routine  each  time  that  it  is  called.  The  counts  are  reset  to  zero 
whenever  a  new  iteration  is  started.  Thus,  if  7  runs  per  iteration 
are  made,  the  acceleration  formula  works  twice;  at  the  end  of  the 
third  run,  and  at  the  end  of  the  sixth  run.  The  subroutine  stores 
the  first  set  of  DFT  at  the  seventh  run.  But  as  at  the  next  ite¬ 
ration,  the  counting  of  the  sets  of  DFT  stored  is  reset  at  0  , 
the  first  set  of  DFT,  which  now  is  that  obtained  at  the  first  run 
of  the  new  iteration,  is  stored  again. 

With  each  test  system  (14,  30,  57  buses),  there  is  a  number  of 
runs  per  iteration  such  that  the  total  time  for  convergence  is  less 
than  the  time  spent  by  Newton-Raphson  method,  provided  that 
a  certain  pair  of  acceleration  factors.  As  an  example  the  following 
values  are  suggested. 

-  for  57  bus  test  system:  6  runs/ite.  along  with  FRL  3.5  and 
FIM  5.0 

-  for  30  bus  test  system:  6  runs/ite.  along  with  FRL  3.0  and 
FIM  5.5 

-  for  14  bus  test  system:  4  runs/ite.  along  with  FRL  3.0  and 
FIM  6.0 

The  results  however  are  liable  to  have  a  small  mismatch  in 
power,  due,  as  when  no  acceleration  factors  were  used,  to  the  cutting 
off  when  the  voltage  increments  are  less  than  the  precision  (0.0001). 

Although  the  solution  is  reached  faster  if  the  suitable  number 
of  runs  per  iteration  and  acceleration  factors  are  used,  the  solution 
still  is  not  good  enough,  because  of  the  slow  power  convergence 
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during  the  final  iterations,  and  the  final  power  mismatch.  In  the 
next  section  it  will  be  shox-m  hox^  a  final  iteration  using  the  Nexvton- 
Raphson  method  reduces  the  time  for  solution  and  eliminates  the  power 
mismatch . 

7.4  ITERATIVE  APPROXIMATE  SOLUTION  AND  SWITCH  TO  AN  EXACT  SOLUTION 

AND  TO  THE  NEWTON-RAPHSON  METHOD 

There  is  no  advantage  in  spending  much  time  on  the  solution  of 
the  linearized  system  during  the  early  iterations  because  it  is  knox\m 
beforehand  that  the  solution  is  only  an  approximation  at  that  stage. 
The  approximate  solution  of  the  linearized  system  in  the  proposed 
method  tries  to  postpone  an  exact,  time-consuming  solution  of  the 
linearized  system  until  the  voltages  are  close  to  the  final  results. 

When  the  results  get  close  enough  to  the  solution,  an  exact  me¬ 
thod  of  solution  can  be  applied  instead  of  the  approximate  iterative 

technique.  The  iteration  in  which  the  change  of  method  of  solu¬ 
tion  of  the  linearized  system  takes  place,  will  be  called  ISW. 

A  procedure  that  greatly  improves  the  performance  of  the  pro¬ 
posed  method  is  the  following:  when  the  iteration  ISW  is  reached, 
the  linearized  system  is  solved  by  triangulation  x^ith  pivoting  instead 
of  solving  it  by  the  Gauss-Seidel  iterative  technique.  From  the 
next  iteration  on,  the  Newton-Raphson  method  is  used.  The  following 
principles  have  to  be  borne  in  mind, 

-  During  the  first  iteration  and  the  early  ones,  only  a  few 
runs  per  iteration  x^ill  be  used.  This  will  give  a  new  set  of  voltages 
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not  as  close  to  the  solution  as  would  have  been  obtained  using 
an  exact  solution.  However,  since  the  time  spent  is  very  short,  the 
process  can  be  repeated  several  times,  and  come  up,  at  the  end,  with 

a  better  approximation  than  that  obtained  with  an  exact  time  consum¬ 
ing  solution.  The  time  spent  for  this  is  less  than  that  for  the 
exact  solution. 

-  The  process  explained  before  is  repeated  until  the  values  of 
the  voltages  are  such  that,  when  applying  an  exact  solution,  it  gives 
a  set  of  values  of  the  voltages  which  go  to  the  final  solution  with 

a  single  iteration  of  the  Newton-Raphson  method.  A  previous  run  of 
the  program  without  any  switching  to  the  Newton-Raphson  method 
could  be  made  to  see  which  is  the  more  convenient  iteration  to  make 
the  switch  to  the  exact  solution.  However  it  is  also  possible  to  let 
the  problem  itself  find  out  when  to  switch  to  the  exact  solution. 

The  maximum  difference  of  power  at  the  beginning  of  the  iteration, 
or  the  maximum  voltage  change  at  the  end  of  it  are  a  good  indication 
of  whether  the  voltages  are  close  or  not  to  the  final  solution. 

-  The  exact  method  used  for  solving  the  linearized  system  at 
the  iteration  ISW  is  triangulation  with  pivoting,  which  means 
destruction  of  the  matrix  of  coefficients.  In  this  case  it  does  not 
matter  because  it  had  to  be  recalculated  with  the  latest  values  of 
the  voltages  so  that  the  Newton-Raphson  technique  can  be  used.  The 
matrix  of  coefficients  dealt  with  was  that  calculated  at  the  first 
iteration  with  the  assumed  voltages,  except  for  the  diagonal  terms, 
which  were  updated  at  each  iteration. 
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-  This  method  of  switching  to  the  Newton-Raphson  technique  is 
one  way  to  improve  the  solution,  but  there  are  several  others.  One, 
that  could  be  tried,  is  increasing  at  each  iteration  the  number  of 
runs  per  iteration.  Both  the  power  mismatch  and  the  slow  convergen¬ 
ce  would  likely  be  reduced  without  a  great  consumption  of  time. 

The  use  of  switching  is  independent  of  acceleration  techniques 
for  solving  the  linearized  system. 

In  figure  7.2  is  displayed;  the  performance  of  the  Newton- 
Raphson  method,  of  the  proposed  method  with  6  runs/ite.  and  FRL  3.5 
and  FIM  5.0  ,  and  of  the  same  method  with  switch  at  the  iteration 
ISW  =  4,  on  the  57  bus  test  system.  In  table  7.4  are  shown  the  times 
and  iterations  spent  by  each  method. 

TABLE  7.4  NEWTON-RAPHSON  METHOD,  PROPOSED  METHOD  WITH  6  RUNS/ITE. 

AND  FRL  3.5  AND  FIM  5.0,  AND  THE  SAME  WITH  SWITCH  AT 
ISW  =  4.  TIMES  AND  ITERATIONS  FOR  SOLUTION  OF  THE  57 
BUS  TEST  SYSTEM 

Solution  by  the  Newton-Raphson  method  using  triangulation  with 
pivoting  for  the  linearized  system.  Total  time  57.8403  s. 

iteration  time  per  new  currents  new  terms  of  sys .  solution 
number  iteration  and  powers  linearized  sys.  and  new  volts. 

1-3  19.2078  0.1592  11.2905  7.7670 


4 


0.1492 


0.1492 


0.0 


0.0 


Figure  7.2  ERROR  IN  THE  REACTIVE  COMPONENT  OF  THE  INCREMENT  IN  VOLTAGE,  DFF ,  VS.  TIME 
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Proposed  method,  4  runs/ite.  FRL  3.5 

FIM  5.0  . 

Total 

time  56.8966 

iteration 

time  per 

new  currents 

new  terms 

of 

sys.  solution 

number 

iteration 

and  powers 

linearized 

sys . 

and  new  volts 

1 

13.7952 

0.1579 

11.6366 

2.0006 

2-21 

2.1493 

0.1488 

0.0113 

1.9892 

Proposed  method,  4  runs 

/ite.  FRL  3.5 

FIM  5.0  . 

ISW  = 

=  4  . 

Total  time 

44.9425  s. 

iteration 

time  per 

new  currents 

new  terms 

of 

sys.  solution 

numb  er 

iteration 

and  powers 

linearized 

sys . 

and  new  volts 

1 

13.4295 

0.1494 

11.2954 

1.9847 

2-3 

2.1488 

0.1481 

0.0112 

1.9896 

4  (ISW) 

7.7244 

0.1476 

0.0112 

7.5657 

5  (N-R) 

19.3418 

0.1496 

11.2930 

7.8992 

6 

0.1493 

0.1493 

0.0 

0.0 

The  same  N.B.  as  in  table  6.3  applies  here. 


The  effectiveness  of  the  exact  solution  of  the  linear  system 
may  be  appreciated.  It  brings  the  voltage  error  almost  to  zero  in 
a  relatively  small  time.  Newt on-Rap hs on ' s  method  needs  then  only 
one  iteration  to  finish  off.  The  time  saved  with  respect  to  the 
Newton-Raphson  method  is  in  this  case  12.9  seconds  (22.3%).  The 
proposed  method  without  switching  to  an  exact  solution,  plus  the 
Newton-Raphson  iteration,  takes  56.9  seconds  to  converge,  i.e.: 
only  saves  0.9  seconds  in  this  case,  with  respect  to  the  plain 
Newton-Raphson  method,  which  takes  57.8  seconds. 


> 

* 

• 

As  an  example  here  are  values  of  the  number  of  runs  per  ite¬ 
ration  and  the  iteration  at  which  to  solve  the  linear  system  in  an 
exact  way,  and  to  sxvitch  to  the  Newton-Raphson  method,  so  that  an 
important  saving  of  time  is  obtained.  For  these  examples,  no 
acceleration  factors  in  the  solution  of  the  linearized  system  were 
used . 


-  for  the  14  bus  test  system 


(3  runs/ite. 
CISW  =  3 


save  16.00% 


,  (2  runs/ite.  1n  ocy 

and  CT„T,  0  save  19.8% 

Clbw  -  d 


-  for  the  30  bus  test  system 

(2  runs/ite.  aqo/ 
CTr,rT  .  save  25.08% 
(ISW  =  4 


,  (3  runs/ite.  0,  A<T/ 

and  CTOrT  o  save  26.0% 

Clbw  =  3 


for  the  57  bus  test  system 
(6  runs/ite. 


;ISW  =  5 


save  18.85% 


,  (12  runs/ite.  cc<7. 

and  CTflT7  0  save  20.55% 

CISW  =  3 


These  are  examples  that  have  just  been  tested.  Other  values 

could  give  better  results.  In  general,  with  the  initial  values  used 

(E°  =  1.0  ,  F°.  =  -0.3),  the  30  bus  test  system  tends  to  yield  more 

k  k 


convergent  results. 


VIII. 


CONCLUSIONS 


8.1  INCONVENIENCES  AND  ADVANTAGES  OF  THE  PROPOSED  METHOD 

A  method  has  been  described  for  solving  the  load-flow  problems. 

It  has  been  tested  with  the  IEEE  14,  30,  and  57  bus  test  systems. 

With  all  three  systems,  the  proposed  method  gives  the  solution  in 
less  time  than  the  Newton-Raphson  method  using  triangulation  with 
pivoting. 

The  proposed  method  is  convergent  for  an  initial  guess  of 
voltages  if  the  particular  guess  yields  convergence  for  the  Newton- 
Raphson  method. 

When  the  proposed  method  is  mentioned  in  this  section,  it  means 
the  proposed  method  as  it  is  employed  in  section  7.4,  that  is  to  say, 
with  an  iterative  approximate  solution  at  the  early  stages,  and  with 
a  switch  to  an  exact  solution  and  the  Newton-Raphson  method  at  the 
end  of  the  process.  There  may  exist  other  ways  to  implement  the 
proposed  method  yielding  more  time  improvement.  Some  alternative  ways 
are  described  in  section  8.3.  In  this  thesis  it  is  not  intended  to 
come  up  with  the  best  implementation  of  the  proposed  method,  but  an 
attempt  has  been  made  to  show  that  there  are  ways  of  implementating 
the  proposed  method  that  yield  the  solution  in  less  time  than  the 
Newton-Raphson  method. 

Since  the  proposed  method  is  a  direct  method,  it  has  all  the 
advantages  and  disadvantages  of  the  direct  methods  over  the  bus  by 
bus  methods,  (see  section  1.2.1),  so  far  as  required  core  and 
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computer  time  taken  by  the  solution  are  concerned. 

The  main  advantage  of  the  proposed  method  over  the  Newton-Raphson 
method  is  that  the  proposed  method  requires  less  time.  The  proposed 
method  needs  three  extra  vectors,  whose  dimensions  are  2(NB  -  S) , 
where  NB  stands  for  the  number  of  buses  and  S  is  the  number  of  slack 
buses.  In  the  57  bus  test  system  these  three  vectors  mean  less  than 
two  per  cent  of  the  total  core.  Moreover,  two  of  the  three  vectors 
are  needed  as  workspace  in  the  acceleration  technique  for  the  itera¬ 
tive  solution  of  the  linearized  system.  However,  a  different  technique 
may  not  demand  the  use  of  the  three  vectors. 

If  several  load-flow  problems  had  to  be  solved  for  the  same 
power  system,  each  having  a  different  set  of  schedules  powers,  the 
proposed  method  could  be  used  in  a  very  efficient  way,  provided  that 
the  same  set  of  initial  voltages  was  used  for  each  load-flow  problem. 
The  fact  of  not  having  to  change  the  non-diagonal  terms  of  the  matrix 
of  coefficients  of  the  linearized  system,  which  takes  more  than  half 
of  the  time  of  the  iteration  if  the  change  takes  place,  means  that  an 
important  reduction  of  time  can  be  achieved. 

The  switch  to  the  exact  solution  in  the  proposed  method  must 
be  done  at  the  suitable  iteration.  It  has  not  been  proven  that  a 
criteria  based  on  the  value  of  the  maximum  power  mismatch,  or  the 
maximum  voltage  change,  or  both,  will  be  suitable  for  every  power 
system.  Therefore,  when  the  proposed  method  is  applied  to  a  load- 
flow  problem  for  the  first  time,  a  risk  exists  of  switching  before 
it  is  convenient,  so  ending  up  running  more  than  one  itteration  with 


. 
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the  Mewton-Raphson  technique.  This  would  mean  in  general,  more  time 
than  the  Newton-Raphson  method  alone. 

If  the  maximum  pox^er  mismatch  or  maximum  voltage  change  criteria 
fails  to  give  the  right  iteration  for  switching  to  the  exact  solution, 
a  trial  and  error  process  could  find  it. 

Another  inconvenience  lies  in  the  fact  that,  when  using  the 
Gauss-Seidel  technique  for  the  solution  of  the  .linearized  system, 
the  problem  is  more  sensitive  to  the  initial  guess  than  it  is  x^hen 
the  Newton-Raphson  method  is  used. 

8.2  REAL  LOAD-FLOW  PROGRAMS  AND  RESEARCH  LOAD-FLOW  PROGRAMS 

The  research  load-flow  programs  handle  a  small  number  of  buses 
in  comparison  with  the  real  load-flow  programs.  Besides,  in  order 
to  squeeze  the  real-load  flox*  program  into  a  non-prohibitive  amount 
of  memory,  many  logical  instructions  are  used  in  the  program  so  that 
only  the  non-zero  elements  of  the  matrix  of  coefficients  of  the 
linearized  system,  and  bus  impedance  or  admittance  matrices  are 
stored  in  the  memory.  This  makes  the  solution  of  the  load-flox^ 
problem  longer,  therefore  a  compromise  must  be  reached  betx^een  the 
amount  of  core  taken,  and  the  time  for  solution.  The  compromise 
chosen  is  different  Tor  different  programs.  Hence  a  fair  com- 
paraison  is  difficult  to  draxxr  between  programs.  Besides  there  are 
other  special  calculations  such  as  voltage  controlled  buses,  acce¬ 
leration  techniques,  net  exchange,  etc.,  which  make  the  comparison 
even  more  difficult.  T.n  the  research  load-flow  programs,  the  tie 
line  interchange,  tap  changing  under  load  transformers,  and  com- 
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plicated  logical  programming  aimed  at  saving  memory,  are  usually 
disregarded . 

In  the  proposed  method,  which  has  been  tested  with  the  IEEE 
14,  30,  and  57  bus  test  system,  some  logical  instructions  have  been 
introduced  so  as  to  be  able  to  deal  with  the  57  bus  test  system 
without  special  permit,  i.e.:  keeping  the  size  of  the  program  under 
200K  (OS  IBM  360  system) ,  which  is  the  limit  core  capacity  for  an 
individual  user  of  the  computer  of  the  University  of  Alberta.  So, 
the  program  works  only  with  the  non-zero  elements  of  the  bus  ad¬ 
mittance  matrix,  which  is  very  sparse.  This  brings  the  program 
employed  here  closer  to  the  real-load  flow  programs.  The  logical 
modifications  introduced  changed  the  amount  of  time  needed  for  the 
solution,  however  the  increasing  or  decreasing  trends  of  the  per¬ 
centages  of  time  taken  by  each  part  of  the  iteration  remained  the 
same. 

Although  the  results  obtained  using  the  proposed  method  are 
faster  than  those  obtained  with  the  Newton-Raphson  method,  it  can 
not  be  concluded  automatically  that  the  proposed  method  is  better 
because  only  the  experience  with  a  real  load-flow  problem  can 
confirm  the  above  mentioned  conclusion.  As  the  Newton-Raphson  and 
the  proposed  method  are  similar,  it  13  hoped  that,  when  adapting 
the  proposed  method  for  a  real  load-flow  problem,  the  same  or 
similar  accelerating  techniques,  and  other  special  features,  can 
be  used  without  difficulty. 
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8.3  SUGGESTIONS  FOR  FURTHER  RESEARCH  ON  THE  PROPOSED  METHOD 

By  taking  advantage  of  not  having  to  calculate  the  non¬ 
diagonal  terms  of  the  matrix  of  coefficients  of  the  linearized 
system,  many  techniques  aimed  at  improving  the  efficiency  of  the 
approximate  solution  of  the  linearized  system  at  each  iteration, 
could  be  implemented. 

Improvement  can  be  brought  about  by  taking  a  variable  number 
of  runs  per  iteration  and  increasing  the  number  of  runs  as  the 
voltages  get  closer  to  the  final  values.  The  criteria  for  increas¬ 
ing  the  number  of  runs  per  iteration  could  be  a  prestablished  one, 
or  one  based  on  the  values  of  the  maximum  mismatch  of  power  or  the 
maximum  change  of  voltage.  This  would  certainly  increase  the 
efficiency  of  the  calculation,  and  hence  it  would  yield  further 
time  savings. 

The  exact  solution  of  the  linearized  system  by  the  triangula¬ 
tion  with  pivoting  technique,  which  is  a  destructive  method,  when 
switching  from  an  approximate  solution  technique  to  an  exact  one, 
could  be  replaced  by  a  non  destructive  technique.  As  a  non-des¬ 
tructive  technique,  the  Gauss-Seidel  iterative  method  could  be  used 
The  exact  solution  may  take  more  time  than  triangulation  with  pivot 
ing,  because  of  the  slow  convergence  of  the  Gauss-Seidel  technique 
for  the  type  of  equation  that  the  linearized  system  is  represented 
by.  On  the  other  hand,  the  coefficients  of  the  linearized  system 
would  not  be  destroyed,  and  two  successive  iterations  with  an  exact 
solution  would  most  likely  lead  to  the  final  solution. 
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The  acceleration  formula  for  the  iterative  solution  of  the 
linearized  system,  can  be  made  more  efficient,  mainly  after  the 
early  iterations,  when  the  results  of  succesive  runs  do  not 
oscillate.  If  the  acceleration  formula  were  more  efficient,  further 
savings  of  time  would  be  obtained. 

Either  triangulation  with  pivoting  or  the  Gauss-Seidel  tech¬ 
nique,  could  be  used  for  the  exact  solution  of  the  linearized  system 
after  the  switch  from  the  approximate  solution.  Then  the  coef¬ 
ficients  of  the  linearized  system  could  be  recalculated,  regarding 
the  last  calculated  voltages  as  the  new  initial  guesses.  Once  this 
is  done,  the  linearized  system  is  solved  with  the  Gauss-Seidel 
technique  using  as  many  iterations  as  necessary  to  guaratee  an 
accuracy  within  the  tolerance. 

The  process  described  above  is  equivalent  to  the  switch  to  an 
exact  solution,  and  an  iteration  of  the  Newton-Raphson  method.  In 
the  suggested  way  however,  at  the  end  of  the  problem,  a  new  matrix 
of  coefficients  of  the  linearized  system  is  available.  This  matrix 
has  been  calculated  with  an  initial  guess  that  is  much  closer  to 
the  solution  than  the  one  that  was  used  before. 
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APPENDIX  1.  AN  ITERATIVE  METHOD  FOR  THE  SOLUTION  OF  A  SET  OF 
QUADRATIC  NONLINEAR  EQUATIONS 

For  the  purpose  of  illustration,  consider  a  very  simple  two 
dimensional  example. 

If  we  have  a  set  of  quadratic  nonlinear  equations  (Al.l) 

P  =  Aef 

Q  =  Bef  (Al.l) 

where  e  and  f  are  the  unknowns,  and  P,Q,A,B,  are  known  quantities, 
a  new  set  of  equivalent  equations,  which  will  simplify  the  applica¬ 
tion  of  an  iterative  way  of  solution,  can  he  derived. 

Let  e*  and  f*  he  the  solution  of  the  system  of  equations 

(Al.l) ,  and  let  e°  and  f°  he  a  pair  of  values  which  are  approxima¬ 

tions  to  the  solutions.  Then  we  can  write 

e*' =  e°  +  Ae° 

f*  =  f°  +  Af°  (A1.2) 

The  solution  of  the  system  (Al.l)  is  e* ,  f*,  which  is  equal 
to  a  starting  value  e°,  f ° ,  plus  the  increments  Ae° ,  Af° ,  which 
brings  the  starting  value  to  the  true  solution  hy  (A1.2). 

If  e*  =  e°  +  Ae°,  and  f*  =  f°  +  Af ° ,  are  substituted  in  (Al.l), 
this  gives, 

P  =  A(e°  +  Ae°)(f°  +  Af°)  =  Ae°f°  +  Ae°Af°  +  Af°Ae°  +  AAe°Af° 

Q  =  B(e°  +  Ae°)(f°  +  Af°)  =  Be°f°  +  Be°Af°  +  Bf°Ae°  +  BAe°Af° 


(A1.3) 
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Calling  P°  =  Ae°f°,  Q°  =  Be°f° ,  and  AP  =  P  -  P° ,  AQ  =  Q  -  Q° , 
and  gathering  terms  in  the  right  hand  side  in  one  of  the  possible 
ways,  we  get  (A1.4) 

P  -  P°  =  AP  =  A(f°  +  Af°) Ae°  +  Ae°Af° 

Q  -  Q°  =  AQ  =  Bf °Ae°  +  B(e°  +  Ae°)Af°  (A1.4) 

where  the  only  unknowns  are  Ae°  and  Af°. 

The  following  iterative  procedure  can  be  used  for  finding  Ae° 
and  Af°.  At  first  we  drop  Ae°  and  Af°  before  e°  and  f°  in  the  right 
hand  side  of  the  equations  (A1.4).  It  results  in  a  linear  system 

(A1.5) 

AP  =  Af °Ae  +  Ae°Af 

AQ  =  Bf °Ae  +  Be°Af  (A1.5) 

Once  the  system  is  solved,  we  obtain  Ael,  Af^,  which  are  not 
the  exact  values  of  the  solution  Ae°,  Af° ,  because  (A1.5)  is  only 
an  approximation  to  the  original  system.  If  we  now  substitute  Ae^, 
Afl  back  in  the  diagonal  coefficients  of  (A1.4)  we  get  a  new  linear 
system  whose  coefficients  are  closer  to  those  of  (A1.4)  than  those 
of  (A1.5),  so  an  improved  solution  can  be  expected  from  the  new 
linearized  system.  The  solution  of  the  new  system,  Ae^,  Af  2 ,  is 
substituted  back  again  in  the  diagonal  coefficients  of  (A1.4),  and 
the  procedure  is  repeated  until  the  improvement  in  the  values  of 


Ae ,  Af,  is  less  than  a  certain  error. 

In  general,  at  the  "ith"  iteration. 


. 
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AP  =  A(f°  +  Afi)  Ae  +  Ae°Af 

AQ  =  Bf°Ae  +  B(e°  +  Ae* )Af  (A1.6) 

and  the  result  would  be  Ae^+1,  Af:*-+~'',  which  would  be  closer  to 
Ae° ,  Af°  than  Ae^  ,  Af i ,  because  the  coefficients  of  Ae  and  Af 
in  ( A1 . 6 )  are  closer  to  their  true  values,  A(f°  +  Af°)  , 

B(e°  +  Ae°)  ,  than  those  used  in  the  former  iteration. 

It  is  obvious  that  if  we  had  (A1.7)  instead  of  (Al.l), 

P  =  Apee  +  A2ef  +  A3ff 

Q  =  Bpee  +  B2ef  +  B3ff  (A1.7) 

or  any  number  of  equations,  with  the  same  number  of  unknowns,  we 
would  always  be  able  to  apply  the  same  method  of  solution,  provided 
that  the  equations  were  quadratic.  See  section  3.1.1  . 


APPENDIX  2.  THE  CALCULATION  OF  THE  Ybus  MATRIX 


For  each  test  power  system,  we  have  as  data,  the  line  impedance, 
and  the  line  charging  admittance.  We  have  to  compute  the  bus  ad¬ 
mittance  matrix,  Y^g . 

There  are  some  transformers  with  off-nominal  turn  ratios ,  and 
synchronous  condensers  connected  to  some  of  the  buses.  The  synchro¬ 
nous  condensers  can  be  replaced  by  a  susceptance. 

The  formulas  are , 
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where  z-^j  are  the  line  impedance,  yij  ,  the  bus  admittances, 
y’ij/2  is  the  one  half  of  total  line  charging  admittances,  and  cf 
the  synchronous  condenser  susceptance,  if  it  exists,  at  the  bus 
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If  between  bus  "p"  and  bus  "q" ,  there  is  a  transformer  with 
an  off-nominal  turn  ratio  "a",  then, 
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The  subroutine  used  for  the  calculations  yields  only  the  non¬ 
zero  elements  of  Y-^g  that  are  above  the  diagonal,  with  a  double 
array  of  numbers  that  give  the  ordered  numbers  of  pairs  of  buses 
that  correspond  to  the  nonzero  elements  of  Ybus* 


APPENDIX  3.  THE  SUBROUTINE  FOR  THE  IMPLEMENTATION  OF  THE  PROPOSED 
METHOD  AND  THE  NEWTON-RAPHSON  METHOD 
Before  the  flow-charts  of  the  Newton-Raphson  method  and  the 
proposed  method,  we  will  give  a  list  of  symbols: 


pok>  Q°k 

4ek»  4pk 

4ek»  Apk 


Ae°k>  Af°k 


AP°k,  AQ°k 


CIRk,  CIIk 
CIR°k ,  CII° 


DCR°k,  DCI°k 


j  ik 


initial  or  estimated  voltage  at  node  Mk" 
powers  calculated  with  the  voltages  e^0,  f-^0 
according  to  the  formula  (1.5) 
increment  of  voltage  between  two  successive 
iterations  of  the  Newton-Raphson  method 
difference  between  the  voltage  at  the  end  of 
an  iteration  of  the  proposed  method  and  the 
initial  voltage  e°k,  f°k 

difference  between  the  true  value  of  the  volt¬ 
ages  and  the  estimated  ones 

difference  between  the  scheduled  power  and  that 
calculated  with  the  initial  values  of  the  voltage 

difference  betweem  the  shheduled  power  and  that 
calculated  with  the  iterated  voltages 
current  at  the  bus  Mk£,  when  the  voltages  are  ek, 
current  at  the  bus  ”kn ,  when  the  voltages  are 

p°  .  f° 

e  1 »  1  i 

difference  between  the  true  value  of  the 
current  at  bus  ”kn  and  its  initial  estimate 
CIR°k,  CII°k 

element  of  the  Jacobian  in  the  Newton-Raphson 


method 
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ik  element  of  the  linearized  system  in  the  proposed  method 

SMX  maximum  mismatch  in  power  with  respect  to  the  scheduled 

powers 

TMX  maximum  change  in  voltage  "between  the  actual  and  the 
former  iteration 
EPS  tolerance  or  precision 

ITE  number  of  the  iteration 

NR  constant  which  is  1  if  the  operations  are  made  as  in  the 

Newton-Raphson  method,  and  0  if  as  in  the  proposed  method 
"true  value"  of  a  variable  stands  for  the  value  that  this  variable 
has  in  the  power  system  when  it  is  operating  with  the  scheduled 
powers . 

Here  follow  some  comments  on  the  flow-charts  and  the  program  used. 

-  A  separate  flow  chart  is  given  for  the  Newton-Raphson  method  and  for 
the  proposed  method  so  that  the  differences  can  be  clearly  seen. 

However  keeping  NR  equal  to  1,  the  flow  chart  for  the  proposed  method 
also  works  for  the  Newton-Raphson  method. 

-  The  voltage  controlled  buses  are  handled  in  the  proposed  method,  in 
the  same  way  as  in  the  Newton-Raphson  method.  Details  can  be  found  in 
section  8.6  of  (3)* 

-  With  the  subroutine  EDTM,  which  is  shown  below,  we  can  effect  either 
the  Newton-Raphson  method  or  the  proposed  method.  Any  bus  can  be  the 
slack  bus,  and  we  can  use  any  number  of  slack  buses.  (See  appendix  6). 

-  To  implement  the  proposed  method  all  we  have  to  do  is  to  keep  NR 


zero  after  the  first  iteration. 
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FLOW  CHART  OF 
THE  METHOD  OF 
NEWTON-RAPHSON 
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SUBROUTINE  ED  iM  (NJ,  NB,  EPS,  FLUX,  LL,  I  N,  JN,  E,  F,  KNS/NZ/  I  OR,  QX,  QN,  J,  G, 
IB, PES/QES/DP/DQ/DCR/DC I ,DFS, DFT/ ERG, ATR, LEL, GRF, LRF, MS, ESR, FSR, 

2 E AC, FAC, ! ORL, JORL, JOR, INF, I NFB, KPT, DPT, EF, FF, DEF, DFF, NGB, IGB, JGB, 

3 NT  I , I CV,DY1,DY2, FAT, FTA) 

REAL*8  J, G, B, PES, QES, DP, DQ, E, F, DCR, DC  I , DFS , DFT, QX, QN, ERG, CP, CM, 

1 E  PS, ESR, FSR, EAC, FAC, PP, PS, AUX, ZER, TMX, SMX 
COMPLEX* 1G  EUX, UEX, ZEC, FLUX, GRF, ATR, PRD 

D I  MENS  I  ON ) J ( NJ, NJ ) , G ( NGB ) , B ( NGB ) , PES ( MB ) , QES ( NB ) , DP ( NB ) ,  DQ( NB ) , 

IE ( NB ) ,  F (NB ) , DCR ( NB ) , DC  I ( NB ) , DFS ( NJ ) , DFT ( NJ ) , QX (NZ ) , QN ( NZ) , 

2ERG ( NZ ) , GRF ( LRF ) , NS ( LRF ) , FLUX ( LL ) , ATR( LEL ) ,  JOR ( NJ ) ,  IOR(NZ), 

3  I  OR  L  (  LED,  JORL  (L  EL),  I  N(LL),  JN(LL),  ESR  (  NB  ) ,  FSR  ( NB  ) ,  EAC  ( NB ) ,  FAC  ( NB )  , 
4 IGB (NGB), JGB (NGB) 

REAL *3  DPT, EF, FF, DEF, DFF, DY1, DY2, FAT, FTA 

D I  MENS  I  ON  DPT ( NJ ) , EF ( NB ) , FF ( NB ) , DEF ( NB ) , DFF ( NB) , DY1 ( NJ ) , DY2 (NJ ) 
DIMENSION  TEMPS ( 100, 3 ) 

350  FORMAT  (6E13.6) 

365  FORMAT  ( //// 13X, 4H I TE  ,12) 

366  FORMAT  ( 12X, I  2, 1H-, I  2, 2X, F10  .  5,  3H  + ( , F10. 5, 2H) J ) 

367  FORMAT  ( /// 11X,  'FINAL  RESULT  7/ 13X, ' BUS ', 12X, 1 P-JQ ',  13X,  '  I  TE= ' , 
112, 2X, 'TIME  USED  * , F8 . 4, '  SECONDS'/) 

368  FORMAT  ( / 13X, 40H I MPOSS I R LE  TO  HOLD  THE  VOLTAGE  IN  NODE  I  2 , 3 X , 
13HQP=,F10. 5/13X, 46HFROM  NOW  ON,  THIS  BUS  WILL  HAVE  CONSTANT  POWER) 

369  FORMAT  ( ///23X, 25H VOLTAGE  CONTROLLED  BUSES  //12X, 3HBUS, 5X, 4HECNT, 
19X, 4HQMAX, 9X, 4HQM IN/) 

370  FORMAT  (12X, I  2, 2X, F10. 5, 3X, F10. 5, 3X,  F10. 5)  I 

371  FORMAT  (//////  13X,  '*****************' //13X,  '  LOAD  FLOW  PROBLEMV/13 
XX, 17 H ******* **********) 

372  FORMAT  (///13X,25H VOLTAGE  CONTROLLED  BUSES  / / 1 2 X , 3HBUS, 5X, 4HECNT, 
111X, 1HQ) 

373  FORMAT  ( 12X, I  2, 2X, F 10 . 5, 3X, F10. 5 ) 

374  FORMAT  (// 13X, 42H***  NO  RESTRICTION  ON  REACTIVE  POWER  ***) 

375  FORMAT  (// 13X, 59H***  RESTRICTIONS  ON  REACTIVE  POWER  TAKEN  INTO  AC 
1COUMT  ***) 

376  FORMAT  ( / 1 IX, 5HS LACK , IX ,  I  2 , F 1 0 . 5 , 3H  + ( , F 10 . 5, 2H ) J ) 

377  FORMAT  ( / 1 IX, 6H LOSSES2X, F 10 . 5, 3H  + ( F 10 . 5, 2H ) J ) 

373  FORMAT  (//13X, 'TIME  USED  ' , F 8 . 4, '  SECONDS') 

379  FORMAT  ( / / 1 3X,  ' BUS ' , 12 X,  '  G- JB  '  , / ) 

330  FORMAT  ( / / 13X, ' SMX  ' , F 10 . 5 ) 

331  FORMAT  (// 13X, ' TMX ' , F10 . 5 ) 

382  FORMAT  ( / // 9X,  ' I TERAT I  ON  ' , IX,  'TIME  PER', 3X, 'NEW  CURRENTS ', 1 X,  ' NEW  i 
1TERMS ' , 3X, 'LINEAR  SYS.  ' , 2  X ,  'TOTAL  TIME'/lOX,  'NUMBER', 3X, ' I TERAT I  Of 
2 ' , 3  X ,  ' AND  POWERS ', IX, ' I N  JACOB  I  AN ', IX, ' AND  NEW  VOLS.  SECONDS') 

383  FORMAT  ( / 12X, I  2 , IX, 5 ( F 12 . 4 ) , IX ) 

Z  E  R  =  0 . 0 

ZEC=(0. 0,0.0) 

KZ  =  0 
KU  =  1 
KD  =  2 
NTT  =  0 


KF  I  =0 

MT=240*( 1000000/ 26) 

WRITE  (6,371) 

WRITE  (6,379) 

WRITF  (6,  366)  (  (  I  GB ( I ) , JGB (  I  ) ,  G  (  I  ) , B ( I ) ) , I = 1, NGB ) 
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55 

IF  ( 1 NF. EQ. 0)  GO  TO  501 

56 

CALL  EDVR  ( KU, NB, NB/ NS, LRF, PES, QES, 4H  BUS,4H 

57 

CALL  EDVR  ( KU, NB, NB, NS, LRF, E, F, 4H  BUS, 4H  EES 

53 

IF  (KNS.EQ.O)  GO  TO  501 

59 

WRITE  (6,369) 

60 

WRITE  (6,370)  ( (  1  0  R  ( l< ) ,  E  R  G  (  K  ) ,  Q  X  (  K  ) ,  Q  N  (  K )  ) ,  K 

61 

WRITE  (6,374) 

62 

501 

1  TE  =  0 

63 

K  D  D  =  0 

64 

IF  (KPT.EQ.l)  KDD= 1 

65 

CALL  CS019A  (MT) 

66 

499 

NR  =  1 

67 

DO  500  1 =1, NB 

68 

EAC ( ! ) =E ( 1 ) 

1  69 

F  A  C ( 1 )=F( 1 ) 

70 

ESR ( 1 ) =E ( ! ) 

71 

500 

FSR ( 1  )=F( 1 ) 

72 

502 

1 TE= 1 TE+ 1 

73 

IF  (INFB.NE.O)  WRITE  (6,365)  1 TE 

74 

IF  (ITE.GE.ICV)  KDD=0 

75 

SMX=Q. 0 

76 

L  1  =  1 

77 

DO  505  1=1, NB 

78 

IF  (1 . NE . NS ( L 1 ) )  GO  TO  503 

79 

IF  (LI.LT.LRF)  LI=LI+1 

80 

GO  TO  506 

81 

503 

PP=0. 0 

32 

PS  =  0. 0 

33 

CP=0. 0 

8  4 

CN=0 . 0 

35 

KSY  =  0 

86 

DO  504  K=l, NGB 

87 

IF  ( 1 GB(K) .GT. 1 .AND. KSY. EQ. 1)  GO  TO  536 

88 

IF  ( 1 . NE. JGB(K) )  GO  TO  537 

89 

PP  =  PP+  E  (  1  G  B  (!<))*  G  (!<)  +  F  (  IGB(K)  )*B(K) 

90 

PS=PS+F ( 1 GB ( K ) )*G(K)-E( IGB(K) )*B(K) 

91 

CP=CP+EAC( IGB(K) )*G(K)+FAC( IGB(K) ) *B ( K ) 

92 

C  N  =  C  N  +  F  A  C ( 1 GB(K) )*G(K)-EAC( IGB(K) )*B(K) 

93 

GO  TO  504 

94 

537 

IF  (  1  .NE. IGB(K) )  GO  TO  5  04 

95 

PP=PP+E( JGB(K) )*G(K)+F( JGB(K) )*B(K) 

96 

PS  =  PS  +  F  ( >JGB  (  K )  )*G(K)-E(JGB(K)  )*B(K) 

97 

CP  =  C  P+ EAC ( JGB ( K ) ) *G ( K ) +FAC ( JGB ( K ) ) *B ( K ) 

98 

CN=CN+FAC(JGB(K) )*G(K) -EAC( JGB (K) )*B(K) 

99 

KS  Y  =  1 

100 

5  04 

CONTINUE 

101 

536 

DCR ( 1 )=CP 

10  2 

DC  1 ( 1 ) =  C  N 

103 

AUX  =  PES( 1  )-E( 1 ) *  PP- F ( 1  )  *  PS 

104 

IF  (NR. EQ. 1)  DP( 1 ) =AUX 

105 

IF  (DABS(AUX) .GT.SMX)  SMX=DABS ( AUX ) 

106 

AUX=QES( 1 ) -F ( 1 ) *PP+E (1  ) *PS 

107 

IF  (KNS.EQ.O)  GO  TO  533 

108 

DO  505  K= 1, NZ 

PES,4H  QES ) 
4H  FES) 
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10  9 

IF  ( 1 .  E  Q . 1 0 R ( K ) )  AUX  =  ERG(K)**2-E( 1 )**2-F( 1 )**2 

110 

IF  (1  .EQ.  1  OR  (  K )  )  DQO  )=AUX 

111 

505 

CONTINUE 

112 

533 

IF  (NR.EQ.l)  DQ( 1 )=AUX 

113 

IF  (DABS(AUX) .GT.SMX)  SMX=DABS ( AUX ) 

114 

5  06 

CONTINUE 

115 

IF  (INFB.NE.O)  WRITE  (6,380)  SMX 

116 

IF  (SMX. GT. EPS)  GO  TO  534 

117 

CALL  CS019B  (MTR ) 

118 

TEMPS ( 1 TE,  3 ) =  F  LOAT (MT-MTR ) *  2 . 6E- 05 

119 

NTT  =  1 

120 

GO  TO  518 

121 

534 

IF  (INFB.NE.O. AND. NR.EQ.l)  CALL  EDVR  ( KZ, NB , NB, NS, LRF ,  DP,  DQ,  '  BUS 

122 

1/,DELP,/'DELQ’) 

123 

CALL  CS019B  (MTR) 

124 

TEMPS ( 1 TE, 1 ) =F  LOAT (MT-MTR) *2 . 6E-05 

125 

1  =0 

126 

LI  =1 

127 

DO  515  1 1 =1, NB 

128 

IF  (II  .  NE.NSUI  )  )  GO  TO  507 

129 

IF  (LI.LT.LRF)  LI =L I + 1 

130 

GO  TO  515 

131 

507 

1  =  1  +  1 

132 

1 NL=I +NB-LRF 

133 

IF  (NR.NE.l)  GO  TO  512 

134 

DFS ( 1 )=DP( 1 1 ) 

135 

DFS ( 1 NL)=DQ( 1  I  ) 

136 

K  =  0 

137 

LK  =  1 

138 

DO  511  KK=1, NB 

139 

IF  (KK.NE.NS(LK) )  GO  TO  508 

140 

IF  (LK.LT.LRF)  LK  =  LK+ 1 

141 

GO  TO  511 

142 

508 

K  =  K+  1 

143 

KNL=K+NB-LRF 

144 

J ( I , K ) =0 . 0 

145 

J( 1 , KN  L ) =0 . 0 

146 

DO  550  N=l, NGB 

14  7 

IF  ( ( 1 GB(N).NE. 1  1  .OR.  JGB(N) .NE.KK) .AND. ( 1 GB ( N ) . NE . KK. OR . JGB ( N  ) . NE 

148 

111))  GO  TO  550 

149 

J( 1 ,K)=E( I  1  )*G(N)-F( 1  1  ) *  B ( N  ) 

150 

J ( 1 , KNL) =E ( 1  1  ) *R ( N  )  +  F ( 1  1  ) *G ( N ) 

151 

GO  TO  549 

152 

550 

CONTINUE 

153 

549 

IF  (KNS.EQ.O)  GO  TO  510 

154 

MS  =  0 

155 

DO  509  M  =  1 , N Z 

156 

IF  (II .NE. IOR(M) )  GO  TO  509 

157 

J( 1 N L , K ) =  0 . 0 

153 

<J(  INL,KNL)=0.0 

159 

MS  =  1 

160 

509 

CONTINUE 

161 

IF  (MS.EQ.l)  GO  TO  511 

162 

510 

J(  1  NL, K ) =J ( 1 , KNL) 
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163 
16  4 

165 

166 
167 
1C  8 
163 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

130 

131 

132 
183 
134 

185 

186 
187 
138 
183 

190 

191 

192 

193 

194 

195 

196 
19  7 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 


J(  I  NL,  K  N  L )  =  -  J  (  I  ,K) 

511  CONTINUE 

512  J( I,  I )=J( I, | )  +  DCR ( I  I  ) 

0(1/1  NL)  =0  (  I  /  I  ND  +  DCI  (II) 

IF  (KNS.EQ.O)  GO  TO  514 
MS  =  0 

DO  513  M  =  1 / N Z 

IF  (II . NE. IOR(M) )  GO  TO  513 
J(  I  NL/  I ) =ES R ( I  I  )  +E  (  I  I  ) 

J( I NL, I NL) =FSR( I  I  )+F( I  I ) 

MS  =  1 

513  CONTINUE 

IF  (MS.EQ.l)  GO  TO  515 

514  J( I N  L / I )=0( I  ML/  I ) -DC  I ( I  I  ) 

J( I NL/  I NL) =J ( I NL/ I ML)+DCR( I  I ) 

515  CONTINUE 

CALL  CS019B  (MTR) 

TEMPS ( I TE/ 2 ) =FLOAT (MT-MTR) *2 . 6E-05 
IF  (KDD.EQ.O)  CALL  RSLT  (J/DFS/DFT/JOR/NJ) 

IF  (KDD.NE.O)  CALL  GSDD  (J/DFS,DFT/DPT/NJ/NS/LRF/NB/NTI/DY1/DY2/ 
1FAT/ FTA/ I TE) 

TMX=0 . 0 
I  =0 
LI  =1 

DO  517  11=1/ NB 

IF  (I  I . N E . NS (LI))  GO  TO  516 

FAC ( I  I  )=0. 0 

EAC (I  I  ) =0 . 0 

IF  (LI.LT.LRF)  LI=LI+1 
GO  TO  517 

516  1=1+1 

I N  L  = I +NB-LRF 

EAC ( I  I ) =ESR ( I  I  ) +DFT ( I ) - E ( I  I  ) 

FAC( I  I )=FSR( II )+DFT( I NL)-F( I  I ) 

IF  ( DABS ( EAC ( I  I  ) ) .GT.TMX)  TMX  =  DABS ( EAC (II)) 

IF  (DABS (FAC (I  I )). GT.TMX)  TMX  =  DABS ( FAC ( II)) 

E ( I  I ) =ES R ( I  I ) +  D F T ( I ) 

F ( II  ) =  F  S  R ( I  I  ) +DFT ( I NL) 

D E F ( I  I  )=EF(I  I  )  —  E ( I  I  ) 

DFF ( I  I ) =FF (I  I ) —  F ( I  I  ) 

517  CONTINUE 

IF  (INFB.NE.O)  WRITE  (6/381)  TMX 

IF  (INFB.NE.O)  CALL  EDVR  (KZ/NB/NB/NS, LRF/DEF/DFF  /  '  BUS'/'  DEF'/' 
1DFF  '  ) 

CALL  CS019B  (MTR) 

TEMPS ( I TE/ 3) =F LOAT (MT-MTR) *2 . 6E-05 

IF  (TMX.LE.EPS)  GO  TO  518 

IF  (INFB.NE.O)  WRITE  (6,378)  TEMPS(ITE,3) 

NR  =  0 

IF  (KDD.EQ.O)  GO  TO  499 
IF  (KDD.NE.O)  GO  TO  502 

518  WRITE  (6,367)  I TE, TEMPS ( I TE, 3 ) 

DO  521  1=1, LL 

UEX=ZEC 
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217 

213 

219 

220 

519 

221 

222 

223 

224 

225 

520 

226 

227 

228 

521 

229 

230 

231 

232 

233 

234 

235 

236 

237 

522 

238 

523 

239 

240 

524 

241 

242 

243 

525 

244 

245 

246 

551 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

553 

257 

552 

258 

259 

260 

261 

262 

263 

264 

265 

266 

267 

263 

269 

270 

DO  519  K=l/ NGB 

I F  (I N( I ) . EQ. I GB (K) . AMD. JN( I ) . EQ. JG8 ( K ) .OR. I N( I ) . EO. JGB(K) . AND. 
1JN(  I  )  .  EQ.  IGB(K)  )  U  E  X  =  - D CM  P  L X  ( G  ( l<  ) ,  -  B  (  K )  ) 

CONTINUE 

EUX=ZEC 

DO  520  K=1/LEL 

IF  (IN(I).EQ.IORL(K) .AND . JN( I ) . EQ. JORL(K) . OR. I N(  I  ) . EQ.  JORL(K)  .AND  . 
1JN( I ) . EQ. I ORL(K) )  EUX=ATR ( K ) 

CONTINUE 

FLUX ( I )  =  ( E ( I M ( I  )  )**2  +  F ( I N( I )  )  **2 ) * (UEX  +  EUX )  +  UEX*DCMPLX( (-ECIN(I))* 
IE ( JN  C I ))-F(IN(l ) )*F(JN( I )))/(E(JN(l ) ) *  F ( I N ( I ) ) - E ( I N ( I ) )  *  F  ( J  N  (  I ) ) ) ) 
WRITE  (6,366)  I N ( I ) , JN ( I ) , F LUX ( I ) 

PRD=ZEC 

DO  524  L  =  l, LRF 
NS L=NS ( L) 

E  U  X  =  Z  E  C 

DO  522  ! =  N S  L , L L 

IF  (ILK!).  LT.  MSL)  GO  TO  522 

IF  (  INC  I  )  .EQ.NSL)  EIJX  =  EUX  +  F  LUX  (  I  ) 

IF  ( INC  I ) .GT.NSL)  GO  TO  523 


PRD  =  PRD+  EUX 
GRF  CL) =EUX 
DO  525  1=1, NB 

PRD  =  PRD +  D CMP LX ( PES ( ! ), -QES (  I  ) ) 

CONTINUE 

WRITE  (6,377)  PRD 

CALL  EDVR  ( KU, NB, MB, NS, LRF, E, F, 4H  BUS, 4HEF I N, 4HFF I N ) 

WRITE  (6,332) 

DO  552  I =1, I TE 
IF  (I.EQ.l)  TMO  =TEM PS ( 1,3) 

IF  (I.GE.2)  TMO  =T EM PS ( I , 3 ) -TEM PS ( I  - 1,  3 ) 

IF  (NTT.EQ.l.AND. I .EQ. ITE)  GO  TO  553 
TM3  =  TEM  PS ( 1 , 3 ) -TEM  PS ( I , 2 ) 

TM2  =TEMPS ( I , 2) -TEMPS ( I ,1) 

IF  Cl  .  GE  .  2  )  TM1=TEMPS  (  I  ,  D-TEMPSC  I  -1,  3) 

IF  ( ! . EQ. 1 )  TM1=TEMPS (1,1) 

I F( I . LT. I TE.OR. NTT. EQ. 0 )WR I TE (6, 383)  I , TMO, TM1, TM2, TM3, TEMPS C I ,3) 


IF  (KNS.EQ.O)  RETURN 
WRITE  (6,372) 

DO  5  27  K  = 1 , N Z 
IF  CIOR(K).EQ.O)  GO  TO  527 
PP=0. 0 
PS=0. 0 
KS  Y  =  0 

DO  526  1=1, NGB 

IF  ( I GB ( I ) . GT. I  OR (K) . AND . KSY. EQ. 1 )  GO  TO  532 
IF  (lOR(K).NE.JGB(l ))  GO  TO  535 
PP  =  PP+E( I G  B ( I  ))*G( I ) +  F ( I G  B ( I  )  )*B( I  ) 

PS=PS+F( I GB ( I ) )*G( I )-E( I GB( I ) )*B( I ) 

GO  TO  526 
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271 

535 

IF  (1  OR  ( K ) . NE. 1 GB ( 1  ) )  GO  TO  526 

111 

PP=  PP+  E ( JGB ( 1 ) )*G( 1 )+F(JGB( ! ) )*B( 1 ) 

273 

PS=PS+F(JGB( 1 ) )*G( 1 )-E(JGB( 1 ) )*B( 1 ) 

274 

K5Y  =  1 

275 

526 

CONTINUE 

276 

532 

DQ(IOR(K))=F( IOR(K))*PP-E(l OR ( K ) ) * PS-QES ( IOR(K) ) 

277 

AUX=DSQRT ( E ( IOR(K))**2+F(IOR(K))**2) 

273 

WRITE  (6,  373  )  1  OR ( K ) , AUX, DQ( 1  OR ( K ) ) 

279 

527 

CONTINUE 

280 

KFI  =0 

281 

DO  531  K  =  1 , N Z 

282 

IF  CIOR(K).EQ.O)  GO  TO  531 

283 

IF  (DQ(  IOR(K)  KGE.QXCK)  )  GO  TO  52  8 

284 

IF  (DQ( lOR(K) ). LE.QN(K) )  GO  TO  529 

285 

GO  TO  531 

286 

528 

QES ( I 0  R ( K ) ) =OES ( 1  OR ( K ) ) +QX ( K ) 

287 

GO  TO  530 

288 

529 

QES ( 1  OR ( K ) ) =QES ( 1  OR ( K ) ) +QN ( K ) 

289 

530 

WRITE  (6,  368  )  I  OR ( K ) , QES ( 1  OR ( K ) ) 

290 

1  OR  (10=0 

291 

KF  1  =1 

292 

531 

CONTINUE 

293 

IF  (KFI.EQ.O)  RETURN 

294 

WRITE  (6,375) 

295 

IF  (INFB.NE.O)  CALL  EDVR  (KU, NB, MB, NS, LRF, EAC,  FAC,  4H 

2  96 

1  FES) 

297 

GO  TO  501 

298 

END 

BUS, 4H 


EES, 4H 
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APPENDIX  k .  OTHER  SUBROUTINES  USED 


A  number  of  subroutines  are  used  to  perform  special  operations 
during  the  execution  or  the  printing  operations  of  the  load-flow 
problem.  They  are  the  following: 

EDVR  which  prints  two  vectors  in  two  columns  side  by  side,  along  with 
the  numbers  of  the  components  of  the  vectors. 

RSLT  which  solves  a  linear  system  by  means  of  triangulation  with 
pivoting  technique 

GSDD  which  solves  a  linear  system  by  means  of  the  Gauss-Seidel 

technique  with  the  modifications  explained  in  section  5*1»2.1  . 
Any  number  of  runs  per  iteration  can  be  used 
FOAC  which  implements  the  acceleration  technique  explained  in 
section  7*1»1  in  the  solution  of  the  linearized  system  by 
means  of  the  subroutine  GSDD 

All  the  programs  have  been  written  in  Fortran  IV,  and  all  the 
variables  and  workspaces  of  all  the  subroutines  have  adjustable  di¬ 
mensions,  so  that  the  dimensions  of  the  arrays  and  matrices,  have  to 
be  written  only  once  in  the  main  program. 


a  ' 


APPENDIX  5-  THE  SOLUTION  OF  THE  LOAD-FLOW  PROBLEMS 


The  final  voltages,  and  scheduled  powers  of  the  l4 ,  30,  and  57 
bus  test  systems  are  shown  next.  The  values  are  those  obtained  after 
the  trial  and  error  process  by  which  it  is  determined  whether  the 
voltage  magnitude  at  the  voltage  controlled  buses  can  be  held  at  the 
specified  values  by  means  of  the  available  reactive  power  at  these 
buses.  In  section  4.1.1  it  is  explained  that  to  measure  the  propertie 
of  the  method,  only  the  first  stage  of  the  trial  and  error  technique 
is  necessary. 

Only  in  the  case  of  the  57  bus  test  system  are  we  able  to  reach 
the  final  solution  at  the  end  of  the  first  trial,  as  far  as  availabi¬ 
lity  of  reactive  power  at  the  voltage  controlled  buses,  is  concerned. 
In  the  l4  and  30  bus  test  systems  we  need  more  than  one  trial  to 
arrive  at  the  solutions  shown  in  this  appendix. 

We  will  call  the  real  power,  P,  the  reactive  power,  Q,  the  real 
voltage,  E,  and  the  reactive  voltage,  F.  The  voltage  controlled 
buses  will  be  signaled,  and  there  will  be  an  indication  whether  the 
scheduled  voltage  magnitude  can  be  held  or  it  is  impossible  to  do  so 
with  the  available  reactive  power. 


RESULTS  OF  THE  l4  BUS  TEST  SYSTEM 


controlled 

bus 

P 

Q 

E 

F 

vol.  magn. 

slack 

1 

2.33499 

-0.43432 

1.06000 

0.00000 

vl.  cont 

.  2 

0.18300 

-0.18978 

1.01658 

-0.083^0 

1.02000 

vl.  cont 

.  3 

-0.94200 

-0.19000 

0.92885 

-0.20713 

impossible 

4 

-0.47800 

0.03900 

0.96111 

-0.17252 

5 

-0.07600 

-0.01600 

0.97284 

-0.14721 

vl.  cont 

.  6 

-0.11200 

-0.01500 

0.97903 

-0.25454 

impossible 

7 

0.00000 

0.00000 

0.97190 

-0.23418 

vl.  cont 

.  8 

0.00000 

0.00000 

0.98208 

-0.23663 

impossible 

9 

-0.29500 

-0.16600 

0.95986 

-0.26304 

10 

-0.09000 

-0.05800 

0.95418 

-0.26484 

11 

-0.03500 

-0.01800 

0.96227 

-0.26124 

■ 


-  Ap.5,  2 


bus 

P 

Q 

E 

F 

controlled 
vol.  magn. 

12 

-0.06100 

-0.01600 

0.95939 

-0.26659 

13 

-0.13500 

-0.05800 

0.95399 

-0.26663 

14 

-0.14900 

-0.05000 

0.93356 

-0.27793 

losses 

0.14499 

-0.18068 

RESULTS  OF 

THE 

30  BUS  TEST 

SYSTEM 

bus 

P 

Q 

E 

P 

controlled 

r 

vol.  magn. 

slack 

1 

2.60989 

-0.16100 

1.06000 

0.00000 

vl.  cont. 

2 

0.18300 

0.37300 

1.03801 

-0.09988 

impossible 

3 

-0.02400 

-0.01200 

1.01057 

-0.14187 

4 

-0.07600 

-0.01600 

0.99731 

-0.16952 

vl.  cont. 

5 

-0.94200 

0.18104 

0.97829 

-O.25IH 

1.01000 

6 

0.00000 

0.00000 

0.99026 

-0.19923 

7 

0.22800 

-0.10900 

0.97063 

-0.22790 

vl.  cont. 

8 

-0.30000 

0.07576 

0.98746 

-0.21218 

1.01000 

9 

0.00000 

0.00000 

1.01742 

-0.26166 

10 

-0.05800 

-0.02000 

1.00388 

-0.28817 

vl.  cont. 

11 

0.00000 

0.16372 

1.04790 

-0.26950 

1.08200 

12 

-0.11200 

-0.07500 

1.01993 

-0.27890 

vl.  cont. 

13 

0.00000 

0.10422 

1.03307 

-0.28249 

1.07100 

Ik 

-0.06200 

-0.01600 

1.00140 

-0.29069 

15 

-0.08200 

-0.02500 

0.99665 

-0.29115 

16 

-0.03500 

-0.01800 

1.00486 

-0.28571 

17 

-0.09000 

-0.05800 

0.99801 

-0.28957 

18 

-0.03200 

-0.00900 

0.98537 

-0.29984 

19 

-0.09500 

-0.03400 

0.98264 

-0.30256 

20 

-0.02200 

-0.00700 

0.98787 

-0.30065 

21 

-0.17500 

-0.11200 

0.98979 

-0.29247 

22 

0.00000 

0.00000 

0. 99040 

-0.29239 

23 

-0.03200 

-0.01600 

0.98429 

-0.29462 

2k 

-0.08700 

-0.06700 

0.97762 

-0.29562 

25 

0.00000 

0.00000 

0.97581 

-0.28703 

26 

-0.03500 

-0.02300 

0.95676 

-0.28906 

27 

0.00000 

0.00000 

0.98414 

-0.27965 

28 

0.00000 

0.00000 

0.98469 

-0.20931 

29 

-0.02400 

-0.00900 

0.95895 

-0.29488 

30 

-0.10600 

-0.01900 

0.94337 

-0.30609 

losses 

0.17589 

-0.92300 

RESULTS  OF 

THE 

57  BUS  TEST 

SYSTEM 

Q 

E 

F 

controlled 

bus 

P 

vol.  magn. 

slack 

1 

k. 2369k 

1.11801 

1.04000 

0.00000 

vl.  cont. 

2 

-0.03000 

-0.88754 

1.00978 

-0.02094 

1.01000 

vl.  cont. 

3 

-0.01000 

-0.18175 

0.97963 

-0.10276 

0.98500 

k 

0.00000 

-0.21000 

0.97275 

-0.12526 

' 

-  Ap. 5  >  3  - 


bus 

P 

Q 

E 

p 

controlled 

JL 

r 

vol.  magn. 

5 

-0.13000 

-0.04000 

0.96566 

-0.14512 

vl.  cont . 

6 

-0.75000 

0.05432 

0.96879 

-0.14781 

0.98000 

7 

0.00000 

0.00000 

0.97556 

-0.13020 

vl.  cont. 

8 

3.00000 

0.40091 

1.00193 

-0.07847 

1.00500 

vl.  cont. 

9 

-1.21000 

-0.77919 

0.96632 

-0.16318 

0.98000 

10 

-0.05000 

-0.02000 

0.96665 

-0.19579 

11 

0.00000 

0.00000 

0.95862 

-0.17239 

vl.  cont. 

12 

-0.67000 

1.04518 

0.99810 

-0.18448 

1.01500 

13 

-0.18000 

-0.02300 

0.96463 

-0.16671 

14 

-0.10500 

-0.05300 

0.95735 

-0.15767 

15 

-0.22000 

-0.05000 

0.98030 

-0.12369 

16 

-0.43000 

-0.03000 

1.00128 

-0.15607 

17 

-0.42000 

-0.08000 

1.01295 

-0.09568 

18 

-0.27200 

-0.09800 

0.97978 

-0.20344 

19 

-0.03300 

-0.00600 

0.94449 

-0.22203 

20 

-0.02300 

-0.01000 

0.93748 

-0.22416 

21 

0.00000 

0.00000 

0.98310 

-0.22575 

22 

0.00000 

0.00000 

0.98454 

-0.22510 

23 

-0.06300 

-0.02100 

0.98290 

-0.22590 

24 

0.00000 

0.00000 

0.97258 

-0.22590 

25 

-0.06300 

-0.03200 

0.93365 

-0.30650 

26 

0.00000 

0.00000 

0.93441 

-0.21544 

27 

-0.09300 

-0.00500 

0.96184 

-0.19594 

28 

-0.04600 

-0.02300 

0.98008 

-0.18133 

29 

-0.17000 

-0.02600 

0.99558 

-0.17147 

30 

-0.03600 

-0.01800 

0.91189 

-0.30902 

31 

-0.05800 

-0.02900 

0.88304 

-0.31070 

32 

-0.01600 

-0.00800 

0.90088 

-0.30166 

33 

-0.03800 

-0.01900 

0.89849 

-0.30156 

34 

0.00000 

0.00000 

0.93027 

-0.23456 

35 

-0.06000 

-0.03000 

0.93806 

-0.23231 

36 

0.00000 

0.00000 

0.94850 

-0.23014 

37 

0.00000 

0.00000 

0.95806 

-0.22912 

38 

-0.14000 

-0.07000 

0.98809 

-0.22338 

39 

0.00000 

0.00000 

0.95587 

-0.22939 

40 

0.00000 

0.00000 

0.94547 

-0.22981 

4l 

-0.06300 

-0.03000 

0.96637 

-0.24233 

42 

-0.07100 

-0.04400 

0.93131 

-0.25887 

43 

-0.02000 

-0.01000 

0.98984 

-0.19879 

44 

-0.12000 

-0.01800 

0.99526 

-0.20901 

45 

0.00000 

0.00000 

1.02257 

-0.16693 

46 

0.00000 

0.00000 

i.o4oo4 

-0.20438 

47 

-0.29700 

-0.11600 

1.00888 

-0.22392 

48 

0.00000 

0.00000 

1.00275 

-0.22440 

49 

-0.18000 

-0.08500 

1.01016 

-0.23209 

50 

-0.21000 

-0.10500 

0.99559 

-0.23745 

51 

-0.18000 

-0.05300 

1.02726 

-0.22837 

52 

-0.04900 

-0.02200 

0.96071 

-0.19542 

53 

-0.20000 

-0.10000 

0.94884 

-0.20607 

54 

-0.04100 

-o.oi4oo 

0.97559 

-0.20221 

55 

-0.06800 

-0.03400 

1.01253 

-0.19318 

56 

-0.07600 

-0.02200 

0.93065 

-0.26804 

57 

-0.06700 

-0.02000 

0.92480 

-0.27545 

losses 

0.27894 

-2.07599 

APPENDIX  6.  OPERATION  WITH  MORE  THAN  ONE  SLACK  BUS 


One  of  the  features  of  the  subroutine  EDTM  is  the  possibility 
of  using  more  than  one  slack  bus,  either  for  the  Newton-Raphson 
method,  or  for  the  proposed  method.  A  slack  bus,  as  defined  in 
section  1.2,  is  one  where  the  two  components  of  the  voltage  are  kept 
constant.  Once  we  know  the  load-flow  solution,  we  can  compute  the 
generated  or  consumed  power  at  this  bus  by  means  of  the  equation 

(1.5). 

The  calculated  powers  at  the  slack  buses  will  balance  the  gene¬ 
rations  and  loads  at  the  other  buses,  and  the  losses  in  the  power 
system. 

If  we  solve  the  load-flow  problem  with  only  one  slack  bus,  and 
then  we  solve  the  same  problem  with  more  than  one  slack  bus,  using 
as  constant  voltages  at  the  slack  buses  the  voltages  obtained  as  so¬ 
lution  of  the  load-flow  with  one  slack  bus,  we  are  bound  to  obtain 
the  same  generations  or  loads  at  the  slack  buses  and  the  same  losses 
in  the  power  system.  Experimentation  proves  that  to  be  true;  only 
slight  differences  are  found  due  to  numerical  errors. 

Each  set  of  different  constant  voltages  for  the  slack  buses 
brings  about  different  power  generations  or  consumptions  at  the  slack 
buses  and  different  losses  in  the  power  system. 

The  more  slack  buses  are  used,  the  more  numerically  stable  is 
the  problem.  Some  trial  and  error  techniques  could  be  tried  on  nu¬ 
merically  instable  load-flow  by  means  of  the  use  of  more  than  one 


slack  bus. 


